Introduction SATPy to ILWISPy

Ensure that you have installed SATPy and eventual other sitepackages required, see also the code field below. The Notebook was tested using Python version 3.9. If you use Python3.6 some SATPy functionality might not be available! The sample data for this notebook is expected in a directory 'Intro_Satpy2ILWISPy' within this notebook folder.

import datetime
import shutil
from IPython.core.display import HTML
import os
import os.path
import sys
import glob
import numpy as np
import ilwis
from satpy.scene import Scene
from satpy.resample import get_area_def
from satpy import find_files_and_readers
import warnings
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

#use the Python function os.getcwd() to retrieve the current folder used and add a sub-folder
work_dir = os.getcwd() + '/Intro_Satpy2ILWISPy'

 #set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)

C:\Mozambique\notebooks\notebooks_RS/Intro_Satpy2ILWISPy

ilwis.version()

'1.0 build 20230510'

# read unzipped file, here extension *.nat
fnames = glob.glob(work_dir+'/*.nat')
print(fnames)

scn = Scene(reader='seviri_l1b_native',  filenames=fnames)

['C:\\Mozambique\\notebooks\\notebooks_RS/Intro_Satpy2ILWISPy\\MSG4-SEVI-MSG15-0100-NA-20230102124243.447000000Z-NA.nat']

scn.available_dataset_names()

['HRV',
 'IR_016',
 'IR_039',
 'IR_087',
 'IR_097',
 'IR_108',
 'IR_120',
 'IR_134',
 'VIS006',
 'VIS008',
 'WV_062',
 'WV_073']

#get list of functions available
dir(scn)

['__class__',
 '__contains__',
 '__delattr__',
 '__delitem__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getitem__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__iter__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__setitem__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_check_known_composites',
 '_compare_areas',
 '_composite_loader',
 '_compute_metadata_from_readers',
 '_copy_datasets_and_wishlist',
 '_create_reader_instances',
 '_datasets',
 '_dependency_tree',
 '_filter_loaded_datasets_from_trunk_nodes',
 '_generate_composite',
 '_generate_composites_from_loaded_datasets',
 '_generate_composites_nodes_from_loaded_datasets',
 '_get_dataarrays_from_identifiers',
 '_get_prereq_datasets',
 '_get_sensor_names',
 '_get_writer_by_ext',
 '_ipython_key_completions_',
 '_load_datasets_by_readers',
 '_read_dataset_nodes_from_storage',
 '_read_datasets_from_storage',
 '_readers',
 '_remove_failed_datasets',
 '_resampled_scene',
 '_resamplers',
 '_slice_area_from_bbox',
 '_slice_data',
 '_slice_datasets',
 '_sort_dataset_nodes_by_reader',
 '_update_dependency_tree',
 '_wishlist',
 'aggregate',
 'all_composite_ids',
 'all_composite_names',
 'all_dataset_ids',
 'all_dataset_names',
 'all_modifier_names',
 'all_same_area',
 'all_same_proj',
 'attrs',
 'available_composite_ids',
 'available_composite_names',
 'available_dataset_ids',
 'available_dataset_names',
 'coarsest_area',
 'copy',
 'crop',
 'end_time',
 'finest_area',
 'generate_possible_composites',
 'get',
 'images',
 'iter_by_area',
 'keys',
 'load',
 'max_area',
 'min_area',
 'missing_datasets',
 'resample',
 'save_dataset',
 'save_datasets',
 'show',
 'slice',
 'start_time',
 'to_geoviews',
 'to_xarray_dataset',
 'unload',
 'values',
 'wishlist']

#note calibration = counts
scn.load(['VIS006'], upper_right_corner='NE', calibration='counts')
scn.keys()

[DataID(name='VIS006', wavelength=WavelengthRange(min=0.56, central=0.635, max=0.71, unit='µm'), resolution=3000.403165817, calibration=<calibration.counts>, modifiers=())]

scn['VIS006']

<xarray.DataArray 'reshape-36ee5efd2d045c8714c9191af5d5a7c4' (y: 3712, x: 3712)>
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
    acq_time  (y) datetime64[ns] NaT NaT NaT NaT NaT NaT ... NaT NaT NaT NaT NaT
    crs       object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",...
  * y         (y) float64 5.569e+06 5.566e+06 ... -5.563e+06 -5.566e+06
  * x         (x) float64 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
Attributes: (12/17)
    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...
    units:                    count
    wavelength:               0.635 µm (0.56-0.71 µm)
    standard_name:            counts
    platform_name:            Meteosat-11
    sensor:                   seviri
    ...                       ...
    name:                     VIS006
    resolution:               3000.403165817
    calibration:              counts
    modifiers:                ()
    _satpy_id:                DataID(name='VIS006', wavelength=WavelengthRang...
    ancillary_variables:      []

print(scn['VIS006'].attrs['start_time'])
print(scn['VIS006'].attrs['end_time'])

2023-01-02 12:30:10.117010
2023-01-02 12:45:09.951615

print(scn)

<xarray.DataArray 'reshape-36ee5efd2d045c8714c9191af5d5a7c4' (y: 3712, x: 3712)>
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
    acq_time  (y) datetime64[ns] NaT NaT NaT NaT NaT NaT ... NaT NaT NaT NaT NaT
    crs       object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",...
  * y         (y) float64 5.569e+06 5.566e+06 ... -5.563e+06 -5.566e+06
  * x         (x) float64 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
Attributes: (12/17)
    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...
    units:                    count
    wavelength:               0.635 µm (0.56-0.71 µm)
    standard_name:            counts
    platform_name:            Meteosat-11
    sensor:                   seviri
    ...                       ...
    name:                     VIS006
    resolution:               3000.403165817
    calibration:              counts
    modifiers:                ()
    _satpy_id:                DataID(name='VIS006', wavelength=WavelengthRang...
    ancillary_variables:      []

d = np.array(list(scn.values()))

np.shape(d)

(1, 3712, 3712)

vis006 = scn['VIS006']
vis006.attrs['area']

Area ID: msg_seviri_fes_3km
Description: MSG SEVIRI Full Earth Scanning service area definition with 3 km resolution
Projection: {'a': '6378169', 'h': '35785831', 'lon_0': '0', 'no_defs': 'None', 'proj': 'geos', 'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 3712
Number of rows: 3712
Area extent: (-5570248.4773, -5567248.0742, 5567248.0742, 5570248.4773)

vis006_meas = d[0]
vis006_meas_int = vis006_meas.astype(int)
print(vis006_meas_int[1600,1600])

115

print(scn.available_composite_names())

['airmass', 'ash', 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection', 'day_microphysics', 'day_microphysics_winter', 'dust', 'fog', 'green_snow', 'hrv_clouds', 'hrv_fog', 'hrv_severe_storms', 'hrv_severe_storms_masked', 'ir108_3d', 'ir_cloud_day', 'ir_overview', 'ir_sandwich', 'natural_color', 'natural_color_nocorr', 'natural_color_raw', 'natural_color_with_night_ir', 'natural_color_with_night_ir_hires', 'natural_enh', 'natural_enh_with_night_ir', 'natural_enh_with_night_ir_hires', 'natural_with_night_fog', 'night_fog', 'night_ir_alpha', 'night_ir_with_background', 'night_ir_with_background_hires', 'night_microphysics', 'overview', 'overview_raw', 'realistic_colors', 'snow', 'vis_sharpened_ir']

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1023.0, 0.1))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(3712,3712,1))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)

rcNew.array2raster(vis006_meas.flatten())
rcNew.pix2value(ilwis.Pixel(1600,1600))

115.0

plt.figure(figsize=(8, 8))
plt.imshow(vis006_meas, interpolation='nearest', cmap='Greys_r')
plt.show()

# sub map creation
rcSelect = ilwis.do('selection',rcNew,"boundingbox(1500 1500, 2499 2499)")
print(rcSelect.size().xsize)
print(rcSelect.size().ysize)
print(rcSelect.size().zsize)
rcSelect.pix2value(ilwis.Pixel(500,500))

1000
1000
1

288.0

def descriptive_statistics(raster_input):
    print('Raster size: ',raster_input.size())
    print()
    print('Map extent: ',raster_input.envelope())
    print()
    coordSys_input = raster_input.coordinateSystem()
    print('Coordinate system: ',coordSys_input.toWKT())
    print()
    print('Proj4: ',coordSys_input.toProj4())
    print()
    datadef_input = raster_input.datadef()
    print('Data type: ', datadef_input.domain())
    stats_input = raster_input.statistics(ilwis.PropertySets.pHISTOGRAM)
    print()
    print('Minimum: ',stats_input[ilwis.PropertySets.pMIN]) 
    print('Maximum: ',stats_input[ilwis.PropertySets.pMAX]) 
    print('Mean: ',stats_input[ilwis.PropertySets.pMEAN]) 
    print('No pixels - also nodata: ',stats_input[ilwis.PropertySets.pCOUNT]) 
    print('No pixels - no nodata: ',stats_input[ilwis.PropertySets.pNETTOCOUNT]) 
    print('Sum all: ',stats_input[ilwis.PropertySets.pSUM]) 

    print()
    x=[a for (a,b) in stats_input.histogram()][:-1] 
    y=[b for (a,b) in stats_input.histogram()][:-1]

    plt.plot(x,y,label='Raster Map values')
    plt.xlabel('Data Range')
    plt.ylabel('Data Frequency')
    plt.title('Raster Histogram')
    plt.legend()

#execute the definition using the selected spectral channel
descriptive_statistics(rcSelect)

Raster size:  Size(1000, 1000, 1)

Map extent:  -1069643.728580 -1930759.437234 1930759.437234 1069643.728580

Coordinate system:  PROJCS["_ANONYMOUS_1004750",GEOCS["_ANONYMOUS_1004750",DATUM["Unknown Datum,ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["GEOS"],,UNIT[meter,1.0]]

Proj4:  +proj=geos +a=6378137 +rf=298.257223563 +h=35785831

Data type:  value

Minimum:  76.0
Maximum:  997.0
Mean:  233.159325
No pixels - also nodata:  1000000.0
No pixels - no nodata:  1000000.0
Sum all:  233159325.0

b1s = ilwis.do('linearstretch',rcSelect, 1)
b1s = ilwis.do('setvaluerange', b1s, 0, 255, 1)

#execute the definition using the selected spectral channel
descriptive_statistics(b1s)

Raster size:  Size(1000, 1000, 1)

Map extent:  -1069643.728580 -1930759.437234 1930759.437234 1069643.728580

Coordinate system:  PROJCS["_ANONYMOUS_1004750",GEOCS["_ANONYMOUS_1004750",DATUM["Unknown Datum,ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["GEOS"],,UNIT[meter,1.0]]

Proj4:  +proj=geos +a=6378137 +rf=298.257223563 +h=35785831

Data type:  value

Minimum:  0.0
Maximum:  255.0
Mean:  55.473107
No pixels - also nodata:  1000000.0
No pixels - no nodata:  1000000.0
Sum all:  55473107.0

b1s.pix2value(ilwis.Pixel(500,500))

76.0

b1s.store('sub_stretched.mpr')

rcNew2 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew) # remove satpy nodata

rcNew2.store('vis006a.mpr')

scn.available_dataset_names()

['HRV',
 'IR_016',
 'IR_039',
 'IR_087',
 'IR_097',
 'IR_108',
 'IR_120',
 'IR_134',
 'VIS006',
 'VIS008',
 'WV_062',
 'WV_073']

scn = Scene(reader='seviri_l1b_native',  filenames=fnames)
scn.load(['VIS006'])
scn.load(['VIS008'])
scn.load(['IR_016'])

r_band = 'IR_016'
g_band = 'VIS008'
b_band = 'VIS006'
scn.load([r_band, g_band, b_band])

from satpy.composites import GenericCompositor
comp = GenericCompositor('my_rgb')
my_rgb = comp((scn[r_band], scn[g_band], scn[b_band]))

from satpy.writers import get_enhanced_image
plt.figure(figsize=(8,8))
img = get_enhanced_image(my_rgb)
# get DataArray out of `XRImage` object
img_data = img.data
img_data.plot.imshow(vmin=0, vmax=1, rgb='bands')

<matplotlib.image.AxesImage at 0x2287f587d90>

scn2 = Scene(reader='seviri_l1b_native',  filenames=fnames)
scn2.load(['VIS006'])
scn2.load(['VIS008'])
scn2.load(['IR_016'])

vis006 = scn2["VIS006"]
vis008 = scn2["VIS008"]
ir_016 = scn2["IR_016"]

vis006_meas = vis006.values#[:-1:]
vis008_meas = vis008.values#[:-1:]
ir_016_meas = ir_016.values#[:-1:]

print(vis006_meas[1600,1600])
print(vis008_meas[1600,1600])
print(ir_016_meas[1600,1600])

vis006_rot = np.rot90(vis006_meas, k=2, axes=(1, 0))
vis008_rot = np.rot90(vis008_meas, k=2, axes=(1, 0))
ir_016_rot = np.rot90(ir_016_meas, k=2, axes=(1, 0))

print(vis006_rot[1600,1600])
print(vis008_rot[1600,1600])
print(ir_016_rot[1600,1600])

data = np.array([vis006_rot, vis008_rot, ir_016_rot]).flatten()

5.3496265
4.144271
3.9428484
6.1138597
22.51347
19.385675

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.001))
rcNew2 = ilwis.RasterCoverage()
rcNew2.setSize(ilwis.Size(3712,3712,3))
rcNew2.setGeoReference(grf)
rcNew2.setDataDef(dfNum)

rcNew2.array2raster(data)

print(rcNew2.pix2value(ilwis.Pixel(1600,1600,0)))
print(rcNew2.pix2value(ilwis.Pixel(1600,1600,1)))
print(rcNew2.pix2value(ilwis.Pixel(1600,1600,2)))

6.1138596534729
22.513469696044922
19.38567543029785

rcNew2 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew2) # remove satpy nodata
rcNew2.store('msg_3_bands.mpl')

print(rcNew2.size())
print(rcNew2.envelope())
rcNew2coordSys = rcNew2.coordinateSystem()
print(rcNew2coordSys.toWKT())
print(rcNew2coordSys.toProj4())

Size(3712, 3712, 3)
-5570248.477300 -5567248.074200 5567248.074200 5570248.477300
PROJCS["_ANONYMOUS_1004811",GEOCS["_ANONYMOUS_1004811",DATUM["Unknown Datum,ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["GEOS"],,UNIT[meter,1.0]]
+proj=geos +a=6378137 +rf=298.257223563 +h=35785831

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.001))
b006 = ilwis.RasterCoverage()
b006.setSize(ilwis.Size(3712,3712,1))
b006.setGeoReference(grf)
b006.setDataDef(dfNum)

b006.array2raster(vis006_rot.flatten())
print(b006.pix2value(ilwis.Pixel(1600,1600)))
print(b006.size())

6.1138596534729
Size(3712, 3712, 1)

b006 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', b006) # remove satpy nodata
b006.store('b006.mpr')

b006_sub = ilwis.do('select',b006, 'boundingbox(1300 250, 1649 499)')

pnew = ilwis.PropertySets.pMIN|ilwis.PropertySets.pMAX|ilwis.PropertySets.pDISTANCE|ilwis.PropertySets.pDELTA|ilwis.PropertySets.pNETTOCOUNT|ilwis.PropertySets.pCOUNT|ilwis.PropertySets.pSUM|ilwis.PropertySets.pMEAN|ilwis.PropertySets.pMEDIAN|ilwis.PropertySets.pPREDOMINANT|ilwis.PropertySets.pSTDEV|ilwis.PropertySets.pHISTOGRAM
new_stat = b006_sub.statistics(pnew)
print("Minimum:", (new_stat.prop(ilwis.PropertySets.pMIN)))
print("Maximum:", (new_stat.prop(ilwis.PropertySets.pMAX)))
print("Mean:", (new_stat.prop(ilwis.PropertySets.pMEAN)))
print("Median:", (new_stat.prop(ilwis.PropertySets.pMEDIAN)))
print("STD:", (new_stat.prop(ilwis.PropertySets.pSTDEV)))
#print("Predominant:", (new_stat.prop(ilwis.PropertySets.pPREDOMINANT))) # check predominant!!!
print("Count:", (new_stat.prop(ilwis.PropertySets.pCOUNT)))

Minimum: 2.2926974296569824
Maximum: 41.077491760253906
Mean: 11.080568362835475
Median: 11.024553825112797
STD: 4.706383701169935
Count: 87500.0

print(b006_sub.size())

Size(350, 250, 1)

b006_sub.store('sub.mpr')

grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=20 40 40 20, gridsize=600 600, cornerofcorners=yes')

b006res = ilwis.do('resample', b006, grf_LL, 'nearestneighbour')
print(b006.envelope())
print(b006res.envelope())
print(b006res.size())
print(b006res.coordinateSystem())
print(b006res.coord2value(ilwis.Coordinate(30, 30))) #querying a location in lon/lat.
print(b006.coord2value(ilwis.Coordinate(2725000, 2888500))) #querying a location in original projection.
print(b006res.pix2value(ilwis.Pixel(100,100))) #querying a location in column and row
print(b006res.coord2value(ilwis.Coordinate(22.78, 5.034))) 

-5570248.477300 -5567248.074200 5567248.074200 5570248.477300
20.000000 20.000000 40.000000 40.000000
Size(600, 600, 1)
LatLon WGS 84
14.998059272766113
17.195228576660156
2.9614005088806152
0.0

b006res.store('b006res.mpr')

Section on use of colour composites

scn4 = Scene(reader='seviri_l1b_native', filenames=fnames)

cir = 'colorized_ir_clouds'
scn4.load([cir], upper_right_corner='NE')
scn4.show(cir)

scn = Scene(reader='seviri_l1b_native', filenames=fnames)
scn.load(['colorized_ir_clouds'], upper_right_corner='NE')
scn['colorized_ir_clouds']

<xarray.DataArray 'getitem-c0b59d411a2fbc7cdfe72b1c7b406998' (y: 3712, x: 3712)>
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
    acq_time  (y) datetime64[ns] NaT NaT NaT NaT NaT NaT ... NaT NaT NaT NaT NaT
    crs       object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",...
  * y         (y) float64 5.569e+06 5.566e+06 ... -5.563e+06 -5.566e+06
  * x         (x) float64 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
Attributes: (12/20)
    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...
    units:                    K
    wavelength:               10.8 µm (9.8-11.8 µm)
    standard_name:            colorized_ir_clouds
    platform_name:            Meteosat-11
    sensor:                   seviri
    ...                       ...
    modifiers:                ()
    _satpy_id:                DataID(name='colorized_ir_clouds', resolution=3...
    ancillary_variables:      []
    optional_datasets:        []
    prerequisites:            [DataQuery(name='IR_108')]
    optional_prerequisites:   []

from satpy.writers import get_enhanced_image

plt.figure(figsize=(8,8))
img = get_enhanced_image(scn['colorized_ir_clouds'])
# get DataArray out of `XRImage` object
img_data = img.data
img_data.plot.imshow(rgb='bands', vmin=0, vmax=1)

<matplotlib.image.AxesImage at 0x2282b06f910>

data1_rot0 = np.rot90(img_data[0], k=4, axes=(1, 0))
data1_rot1 = np.rot90(img_data[1], k=4, axes=(1, 0))
data1_rot2 = np.rot90(img_data[2], k=4, axes=(1, 0))

data1_rot = np.array([data1_rot2, data1_rot1, data1_rot0]).flatten()

data1 = np.array(data1_rot).flatten()

print(data1.ndim)
print(data1.shape)
print(data1.size)

1
(41336832,)
41336832

print(len(data1))

41336832

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.001))
rcNew3 = ilwis.RasterCoverage()
rcNew3.setSize(ilwis.Size(3712,3712,3))
rcNew3.setGeoReference(grf)
rcNew3.setDataDef(dfNum)

rcNew3.array2raster(data1)
rcNew3 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew3) # remove satpy nodata

print(rcNew3.pix2value(ilwis.Pixel(1500,2400,0)))
print(rcNew3.pix2value(ilwis.Pixel(1500,2400,1)))
print(rcNew3.pix2value(ilwis.Pixel(1500,2400,2)))

0.19432323469837695
0.19439588613339462
0.19442534111468598

print(rcNew3.size())
print(rcNew3.envelope())
rcNew3coordSys = rcNew3.coordinateSystem()
print(rcNew3coordSys.toWKT())
print(rcNew3coordSys.toProj4())

Size(3712, 3712, 3)
-5570248.477300 -5567248.074200 5567248.074200 5570248.477300
PROJCS["_ANONYMOUS_1004927",GEOCS["_ANONYMOUS_1004927",DATUM["Unknown Datum,ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["GEOS"],,UNIT[meter,1.0]]
+proj=geos +a=6378137 +rf=298.257223563 +h=35785831

rcNew3.store('cir.mpl')

Section on use of NWC-SAF data (format netCDF)

myfiles = glob.glob(work_dir+'/S_NWC_CRR_MSG4_global-VISIR_20230117T051500Z.nc')

crr = Scene(filenames=myfiles, reader='nwcsaf-geo')
print(crr.available_composite_names())
print(crr.available_dataset_names())

['convective_precipitation_hourly_accumulation', 'convective_rain_rate']
['crr', 'crr_accum', 'crr_accum_pal', 'crr_conditions', 'crr_intensity', 'crr_intensity_pal', 'crr_pal', 'crr_quality', 'crr_status_flag']

#get listing of other functions
dir(crr)

['__class__',
 '__contains__',
 '__delattr__',
 '__delitem__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getitem__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__iter__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__setitem__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 '_check_known_composites',
 '_compare_areas',
 '_composite_loader',
 '_compute_metadata_from_readers',
 '_copy_datasets_and_wishlist',
 '_create_reader_instances',
 '_datasets',
 '_dependency_tree',
 '_filter_loaded_datasets_from_trunk_nodes',
 '_generate_composite',
 '_generate_composites_from_loaded_datasets',
 '_generate_composites_nodes_from_loaded_datasets',
 '_get_dataarrays_from_identifiers',
 '_get_prereq_datasets',
 '_get_sensor_names',
 '_get_writer_by_ext',
 '_ipython_key_completions_',
 '_load_datasets_by_readers',
 '_read_dataset_nodes_from_storage',
 '_read_datasets_from_storage',
 '_readers',
 '_remove_failed_datasets',
 '_resampled_scene',
 '_resamplers',
 '_slice_area_from_bbox',
 '_slice_data',
 '_slice_datasets',
 '_sort_dataset_nodes_by_reader',
 '_update_dependency_tree',
 '_wishlist',
 'aggregate',
 'all_composite_ids',
 'all_composite_names',
 'all_dataset_ids',
 'all_dataset_names',
 'all_modifier_names',
 'all_same_area',
 'all_same_proj',
 'attrs',
 'available_composite_ids',
 'available_composite_names',
 'available_dataset_ids',
 'available_dataset_names',
 'coarsest_area',
 'copy',
 'crop',
 'end_time',
 'finest_area',
 'generate_possible_composites',
 'get',
 'images',
 'iter_by_area',
 'keys',
 'load',
 'max_area',
 'min_area',
 'missing_datasets',
 'resample',
 'save_dataset',
 'save_datasets',
 'show',
 'slice',
 'start_time',
 'to_geoviews',
 'to_xarray_dataset',
 'unload',
 'values',
 'wishlist']

crr.load(['convective_precipitation_hourly_accumulation', 'convective_rain_rate'])

#crr.load(['convective_precipitation_hourly_accumulation'])

print(crr)

<xarray.DataArray 'where-e3ff70549cdc1eec3dcda261e978624f' (bands: 3, y: 3712, x: 3712)>
dask.array<where, shape=(3, 3712, 3712), dtype=float64, chunksize=(1, 3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
  * bands    (bands) <U1 'R' 'G' 'B'
Attributes: (12/23)
    _satpy_id:               DataID(name='convective_rain_rate', resolution=3...
    resolution:              3000
    _FillValue:              nan
    reader:                  nwcsaf-geo
    coordinates:             lon lat
    flag_meanings:           00_02_mm_h 02_1_mm_h 1_2_mm_h 2_3_mm_h 3_5_mm_h ...
    ...                      ...
    units:                   1
    modifiers:               ()
    standard_name:           convective_rain_rate
    prerequisites:           ['crr', 'crr_pal']
    optional_prerequisites:  []
    mode:                    RGB
<xarray.DataArray 'where-024f6f2f18a9ee52fe51458f101b06d5' (bands: 3, y: 3712, x: 3712)>
dask.array<where, shape=(3, 3712, 3712), dtype=float64, chunksize=(1, 3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
  * bands    (bands) <U1 'R' 'G' 'B'
Attributes: (12/21)
    _satpy_id:               DataID(name='convective_precipitation_hourly_acc...
    resolution:              3000
    _FillValue:              nan
    reader:                  nwcsaf-geo
    coordinates:             lon lat
    area:                    Area ID: some_area_name\nDescription: On-the-fly...
    ...                      ...
    units:                   mm
    modifiers:               ()
    standard_name:           convective_precipitation_hourly_accumulation
    prerequisites:           ['crr_accum', 'crr_pal']
    optional_prerequisites:  []
    mode:                    RGB

crr.show('convective_rain_rate')

crr_load = crr['convective_rain_rate']

crr_meas = crr_load.values

arr_reshaped = np.transpose(crr_meas, (1, 2, 0))

plt.imshow(arr_reshaped, interpolation='nearest')
plt.show()

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.00))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(3712,3712,3))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)

rcNew.array2raster(crr_meas.flatten())
rcNew = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew) # remove satpy nodata
print(rcNew.coord2value(ilwis.Coordinate(4896095.17,-1425952.85))) 

0.7019607843137254

rcNew.store('crr_class.mpr')

crr.load(['crr_accum'])

print(crr)

<xarray.DataArray 'where-e3ff70549cdc1eec3dcda261e978624f' (bands: 3, y: 3712, x: 3712)>
dask.array<where, shape=(3, 3712, 3712), dtype=float64, chunksize=(1, 3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
  * bands    (bands) <U1 'R' 'G' 'B'
Attributes: (12/23)
    _satpy_id:               DataID(name='convective_rain_rate', resolution=3...
    resolution:              3000
    _FillValue:              nan
    reader:                  nwcsaf-geo
    coordinates:             lon lat
    flag_meanings:           00_02_mm_h 02_1_mm_h 1_2_mm_h 2_3_mm_h 3_5_mm_h ...
    ...                      ...
    units:                   1
    modifiers:               ()
    standard_name:           convective_rain_rate
    prerequisites:           ['crr', 'crr_pal']
    optional_prerequisites:  []
    mode:                    RGB
<xarray.DataArray 'where-024f6f2f18a9ee52fe51458f101b06d5' (bands: 3, y: 3712, x: 3712)>
dask.array<where, shape=(3, 3712, 3712), dtype=float64, chunksize=(1, 3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
  * bands    (bands) <U1 'R' 'G' 'B'
Attributes: (12/21)
    _satpy_id:               DataID(name='convective_precipitation_hourly_acc...
    resolution:              3000
    _FillValue:              nan
    reader:                  nwcsaf-geo
    coordinates:             lon lat
    area:                    Area ID: some_area_name\nDescription: On-the-fly...
    ...                      ...
    units:                   mm
    modifiers:               ()
    standard_name:           convective_precipitation_hourly_accumulation
    prerequisites:           ['crr_accum', 'crr_pal']
    optional_prerequisites:  []
    mode:                    RGB
<xarray.DataArray 'crr_accum' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=float64, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
Attributes: (12/16)
    long_name:            NWC GEO CRR Convective Hourly Rainfall Accumulation
    units:                mm
    valid_range:          [ 0. 50.]
    _FillValue:           nan
    ancillary_variables:  [<xarray.DataArray 'crr_status_flag' (y: 3712, x: 3...
    coordinates:          lon lat
    ...                   ...
    reader:               nwcsaf-geo
    area:                 Area ID: some_area_name\nDescription: On-the-fly ar...
    name:                 crr_accum
    resolution:           3000
    modifiers:            ()
    _satpy_id:            DataID(name='crr_accum', resolution=3000, modifiers...

crr.show('crr_accum')

crr_meas = crr['crr_accum'].values

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.000))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(3712,3712,1))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)

rcNew.array2raster(crr_meas.flatten())
rcNew = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew) # remove satpy nodata
print(rcNew.coord2value(ilwis.Coordinate(4896095.17,-1425952.85))) 

22.40000033378601

rcNew.store('crr_value.mpr')

myfiles1 = glob.glob(work_dir+'/S_NWC_RDT-CW_MSG4_global-VISIR_20230117T034500Z.nc')

rdt = Scene(filenames=myfiles1, reader='nwcsaf-geo')
print(rdt.available_composite_names())
print(rdt.available_dataset_names())

['rdt_cell_type']
['MapCellCatType', 'MapCellCatType_pal', 'MapCell_conditions', 'MapCell_quality']

rdt.load(['MapCellCatType', 'MapCellCatType_pal'])

rdt.load(['rdt_cell_type'])

print(rdt)

<xarray.DataArray 'MapCellCatType' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=uint8, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
Attributes: (12/19)
    long_name:            NWC GEO RDT-CW Type and phase of significant cells
    valid_range:          [ 0 21]
    _FillValue:           255
    ancillary_variables:  [<xarray.DataArray 'MapCellCatType_pal' (pal01_colo...
    coordinates:          lon lat                  
    comment:              0:Non_convective 1:Convective_triggering 2:Convecti...
    ...                   ...
    reader:               nwcsaf-geo
    area:                 Area ID: some_area_name\nDescription: On-the-fly ar...
    name:                 MapCellCatType
    resolution:           3000
    modifiers:            ()
    _satpy_id:            DataID(name='MapCellCatType', resolution=3000, modi...
<xarray.DataArray 'MapCellCatType_pal' (pal01_colors: 256, pal_RGB: 3)>
dask.array<add, shape=(256, 3), dtype=uint8, chunksize=(256, 3), chunktype=numpy.ndarray>
Dimensions without coordinates: pal01_colors, pal_RGB
Attributes: (12/15)
    long_name:            RGB palette for MapCellCatType
    valid_range:          [  0 255]
    colormodel:           RGB
    comment:              Palette applicable to field MapCellCatType
    platform_name:        Meteosat-11
    sensor:               {'seviri'}
    ...                   ...
    end_time:             2023-01-17 03:57:22
    reader:               nwcsaf-geo
    name:                 MapCellCatType_pal
    resolution:           3000
    modifiers:            ()
    _satpy_id:            DataID(name='MapCellCatType_pal', resolution=3000, ...
<xarray.DataArray 'where-b77cc868e631873eeb00904cb6ff3808' (bands: 3, y: 3712, x: 3712)>
dask.array<where, shape=(3, 3712, 3712), dtype=float64, chunksize=(1, 3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float64 5.569e+06 5.566e+06 5.563e+06 ... -5.563e+06 -5.566e+06
  * x        (x) float32 -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
    crs      object PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",E...
  * bands    (bands) <U1 'R' 'G' 'B'
Attributes: (12/24)
    _FillValue:              nan
    reader:                  nwcsaf-geo
    name:                    rdt_cell_type
    end_time:                2023-01-17 03:57:22
    platform_name:           Meteosat-11
    _satpy_id:               DataID(name='rdt_cell_type', resolution=3000)
    ...                      ...
    units:                   1
    modifiers:               ()
    standard_name:           rdt_cell_type
    prerequisites:           ['MapCellCatType', 'MapCellCatType_pal']
    optional_prerequisites:  []
    mode:                    RGB

rdt.show('rdt_cell_type')

rdt_meas = rdt['rdt_cell_type'].values

grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=geos +h=35785831 +a=6378137 +rf=298.257223563,envelope=-5570248.4773 5570248.4773 5567248.0742 -5567248.0742,gridsize=3712 3712,cornerofcorners=yes')
dfNum = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0.0, 1000.0, 0.000))
rcNew = ilwis.RasterCoverage()
rcNew.setSize(ilwis.Size(3712,3712,3))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)

rcNew.array2raster(rdt_meas.flatten())

rcNew.store('rdt_class.mpr')
#note to display in ilwis desktop the layers b0, b1 and b2 as RGB