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.

In [1]:
#load required required libraries and install sitepackages
from ctypes import *
lib1 = cdll.LoadLibrary(r'/home/jovyan/.local/lib/libGLdispatch.so.0.0.0')
lib1 = cdll.LoadLibrary(r'/home/jovyan/.local/lib/libGLX.so.0.0.0')
lib1 = cdll.LoadLibrary(r'/home/jovyan/.local/lib/libGL.so.1.7.0')

lib1 = cdll.LoadLibrary('/opt/conda/lib/libpython3.11.so.1.0')
lib1 = cdll.LoadLibrary('/opt/conda/envs/openeo/lib/libQt5Gui.so.5')
lib1 = cdll.LoadLibrary('/opt/conda/envs/openeo/lib/libQt5Sql.so.5')
lib1 = cdll.LoadLibrary('/opt/conda/envs/openeo/lib/libQt5Concurrent.so.5')

!pip install satpy

Requirement already satisfied: satpy in /opt/conda/lib/python3.11/site-packages (0.43.0)
Requirement already satisfied: numpy>=1.13 in /opt/conda/lib/python3.11/site-packages (from satpy) (1.25.1)
Requirement already satisfied: pillow in /opt/conda/lib/python3.11/site-packages (from satpy) (10.0.0)
Requirement already satisfied: pyresample>=1.24.0 in /opt/conda/lib/python3.11/site-packages (from satpy) (1.27.1)
Requirement already satisfied: trollsift in /opt/conda/lib/python3.11/site-packages (from satpy) (0.5.0)
Requirement already satisfied: trollimage>=1.20 in /opt/conda/lib/python3.11/site-packages (from satpy) (1.20.1)
Requirement already satisfied: pykdtree in /opt/conda/lib/python3.11/site-packages (from satpy) (1.3.7.post0)
Requirement already satisfied: pyyaml>=5.1 in /opt/conda/lib/python3.11/site-packages (from satpy) (6.0)
Requirement already satisfied: xarray!=0.13.0,>=0.10.1 in /opt/conda/lib/python3.11/site-packages (from satpy) (2023.7.0)
Requirement already satisfied: dask[array]>=0.17.1 in /opt/conda/lib/python3.11/site-packages (from satpy) (2023.7.0)
Requirement already satisfied: pyproj>=2.2 in /opt/conda/lib/python3.11/site-packages (from satpy) (3.6.0)
Requirement already satisfied: zarr in /opt/conda/lib/python3.11/site-packages (from satpy) (2.15.0)
Requirement already satisfied: donfig in /opt/conda/lib/python3.11/site-packages (from satpy) (0.8.1)
Requirement already satisfied: appdirs in /opt/conda/lib/python3.11/site-packages (from satpy) (1.4.4)
Requirement already satisfied: packaging in /opt/conda/lib/python3.11/site-packages (from satpy) (23.1)
Requirement already satisfied: pooch in /opt/conda/lib/python3.11/site-packages (from satpy) (1.7.0)
Requirement already satisfied: pyorbital in /opt/conda/lib/python3.11/site-packages (from satpy) (1.8.0)
Requirement already satisfied: click>=8.0 in /opt/conda/lib/python3.11/site-packages (from dask[array]>=0.17.1->satpy) (8.1.5)
Requirement already satisfied: cloudpickle>=1.5.0 in /opt/conda/lib/python3.11/site-packages (from dask[array]>=0.17.1->satpy) (2.2.1)
Requirement already satisfied: fsspec>=2021.09.0 in /opt/conda/lib/python3.11/site-packages (from dask[array]>=0.17.1->satpy) (2023.6.0)
Requirement already satisfied: partd>=1.2.0 in /opt/conda/lib/python3.11/site-packages (from dask[array]>=0.17.1->satpy) (1.4.0)
Requirement already satisfied: toolz>=0.10.0 in /opt/conda/lib/python3.11/site-packages (from dask[array]>=0.17.1->satpy) (0.12.0)
Requirement already satisfied: importlib-metadata>=4.13.0 in /opt/conda/lib/python3.11/site-packages (from dask[array]>=0.17.1->satpy) (6.8.0)
Requirement already satisfied: certifi in /opt/conda/lib/python3.11/site-packages (from pyproj>=2.2->satpy) (2023.5.7)
Requirement already satisfied: setuptools>=3.2 in /opt/conda/lib/python3.11/site-packages (from pyresample>=1.24.0->satpy) (68.0.0)
Requirement already satisfied: configobj in /opt/conda/lib/python3.11/site-packages (from pyresample>=1.24.0->satpy) (5.0.8)
Requirement already satisfied: shapely in /opt/conda/lib/python3.11/site-packages (from pyresample>=1.24.0->satpy) (2.0.1)
Requirement already satisfied: pandas>=1.4 in /opt/conda/lib/python3.11/site-packages (from xarray!=0.13.0,>=0.10.1->satpy) (2.0.3)
Requirement already satisfied: platformdirs>=2.5.0 in /opt/conda/lib/python3.11/site-packages (from pooch->satpy) (3.8.1)
Requirement already satisfied: requests>=2.19.0 in /opt/conda/lib/python3.11/site-packages (from pooch->satpy) (2.31.0)
Requirement already satisfied: scipy in /opt/conda/lib/python3.11/site-packages (from pyorbital->satpy) (1.11.1)
Requirement already satisfied: asciitree in /opt/conda/lib/python3.11/site-packages (from zarr->satpy) (0.3.3)
Requirement already satisfied: fasteners in /opt/conda/lib/python3.11/site-packages (from zarr->satpy) (0.18)
Requirement already satisfied: numcodecs>=0.10.0 in /opt/conda/lib/python3.11/site-packages (from zarr->satpy) (0.11.0)
Requirement already satisfied: zipp>=0.5 in /opt/conda/lib/python3.11/site-packages (from importlib-metadata>=4.13.0->dask[array]>=0.17.1->satpy) (3.16.0)
Requirement already satisfied: entrypoints in /opt/conda/lib/python3.11/site-packages (from numcodecs>=0.10.0->zarr->satpy) (0.4)
Requirement already satisfied: python-dateutil>=2.8.2 in /opt/conda/lib/python3.11/site-packages (from pandas>=1.4->xarray!=0.13.0,>=0.10.1->satpy) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.11/site-packages (from pandas>=1.4->xarray!=0.13.0,>=0.10.1->satpy) (2023.3)
Requirement already satisfied: tzdata>=2022.1 in /opt/conda/lib/python3.11/site-packages (from pandas>=1.4->xarray!=0.13.0,>=0.10.1->satpy) (2023.3)
Requirement already satisfied: locket in /opt/conda/lib/python3.11/site-packages (from partd>=1.2.0->dask[array]>=0.17.1->satpy) (1.0.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/conda/lib/python3.11/site-packages (from requests>=2.19.0->pooch->satpy) (3.2.0)
Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.11/site-packages (from requests>=2.19.0->pooch->satpy) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/conda/lib/python3.11/site-packages (from requests>=2.19.0->pooch->satpy) (2.0.3)
Requirement already satisfied: six in /opt/conda/lib/python3.11/site-packages (from configobj->pyresample>=1.24.0->satpy) (1.16.0)
In [2]:
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')

gdal connector: ilwis3 connector: Organizing data:
In [3]:
#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)

gdal connector: /home/jovyan/mystorage/test/notebooks/ILWISPy/Intro_Satpy2ILWISPy
ilwis3 connector: Organizing data:
In [4]:
ilwis.version()
Out[4]:
'1.0 build 20230718'
In [5]:
# read unzipped file, here extension *.nat
fnames = glob.glob(work_dir+'/*.nat')
print(fnames)

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

['/home/jovyan/mystorage/test/notebooks/ILWISPy/Intro_Satpy2ILWISPy/MSG4-SEVI-MSG15-0100-NA-20230102124243.447000000Z-NA.nat']
In [6]:
scn.available_dataset_names()
Out[6]:
['HRV',
 'IR_016',
 'IR_039',
 'IR_087',
 'IR_097',
 'IR_108',
 'IR_120',
 'IR_134',
 'VIS006',
 'VIS008',
 'WV_062',
 'WV_073']
In [7]:
#get list of functions available
dir(scn)
Out[7]:
['__class__',
 '__contains__',
 '__delattr__',
 '__delitem__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getitem__',
 '__getstate__',
 '__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_area_defs',
 '_compare_areas',
 '_compare_swath_defs',
 '_contained_sensor_names',
 '_copy_datasets_and_wishlist',
 '_create_reader_instances',
 '_datasets',
 '_dependency_tree',
 '_filter_loaded_datasets_from_trunk_nodes',
 '_gather_all_areas',
 '_generate_composite',
 '_generate_composites_from_loaded_datasets',
 '_generate_composites_nodes_from_loaded_datasets',
 '_get_finalized_destination_area',
 '_get_prereq_datasets',
 '_get_writer_by_ext',
 '_ipython_key_completions_',
 '_load_datasets_by_readers',
 '_prepare_resampler',
 '_read_dataset_nodes_from_storage',
 '_read_datasets_from_storage',
 '_reader_times',
 '_readers',
 '_reduce_data',
 '_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',
 'chunk',
 'coarsest_area',
 'compute',
 'copy',
 'crop',
 'end_time',
 'finest_area',
 'generate_possible_composites',
 'get',
 'images',
 'iter_by_area',
 'keys',
 'load',
 'max_area',
 'min_area',
 'missing_datasets',
 'persist',
 'resample',
 'save_dataset',
 'save_datasets',
 'sensor_names',
 'show',
 'slice',
 'start_time',
 'to_geoviews',
 'to_xarray',
 'to_xarray_dataset',
 'unload',
 'values',
 'wishlist']
In [8]:
#note calibration = counts
scn.load(['VIS006'], upper_right_corner='NE', calibration='counts')
scn.keys()
Out[8]:
[DataID(name='VIS006', wavelength=WavelengthRange(min=0.56, central=0.635, max=0.71, unit='µm'), resolution=3000.403165817, calibration=<4>, modifiers=())]
In [9]:
scn['VIS006']
Out[9]:
<xarray.DataArray 'reshape-e4132bd22c06954dad3e0aad9b837671' (y: 3712, x: 3712)>
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(928, 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/18)
    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...
    time_parameters:          {'nominal_start_time': datetime.datetime(2023, ...
    units:                    count
    wavelength:               0.635 µm (0.56-0.71 µm)
    standard_name:            counts
    platform_name:            Meteosat-11
    ...                       ...
    name:                     VIS006
    resolution:               3000.403165817
    calibration:              counts
    modifiers:                ()
    _satpy_id:                DataID(name='VIS006', wavelength=WavelengthRang...
    ancillary_variables:      []
xarray.DataArray
'reshape-e4132bd22c06954dad3e0aad9b837671'
  • y: 3712
  • x: 3712
  • dask.array<chunksize=(928, 3712), meta=np.ndarray>
    Array Chunk
    Bytes 52.56 MiB 13.14 MiB
    Shape (3712, 3712) (928, 3712)
    Dask graph 4 chunks in 35 graph layers
    Data type float32 numpy.ndarray
    		3712 3712
    • acq_time
      (y)
      datetime64[ns]
      NaT NaT NaT NaT ... NaT NaT NaT NaT
      long_name :
      Mean scanline acquisition time
      array(['NaT', 'NaT', 'NaT', ..., 'NaT', 'NaT', 'NaT'],
      dtype='datetime64[ns]')
    • crs
      ()
      object
      PROJCRS["unknown",BASEGEOGCRS["u...
      array(<Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unk ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Geostationary Satellite (Sweep Y)
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich
, dtype=object)
    • y
      (y)
      float64
      5.569e+06 5.566e+06 ... -5.566e+06
      units :
      meter
      array([ 5568748.275757,  5565747.872591,  5562747.469425, ..., -5559747.066259,
       -5562747.469425, -5565747.872591])
    • x
      (x)
      float64
      -5.569e+06 -5.566e+06 ... 5.566e+06
      units :
      meter
      array([-5568748.275757, -5565747.872591, -5562747.469425, ...,  5559747.066259,
        5562747.469425,  5565747.872591])
    • y
      PandasIndex
      PandasIndex(Index([ 5568748.275756836,  5565747.872591019,  5562747.469425201,
        5559747.066259384,  5556746.663093567,   5553746.25992775,
        5550745.856761932,  5547745.453596115,  5544745.050430298,
        5541744.647264481,
       ...
       -5538744.244098663, -5541744.647264481, -5544745.050430298,
       -5547745.453596115, -5550745.856761932,  -5553746.25992775,
       -5556746.663093567, -5559747.066259384, -5562747.469425201,
       -5565747.872591019],
      dtype='float64', name='y', length=3712))
    • x
      PandasIndex
      PandasIndex(Index([-5568748.275756836, -5565747.872591019, -5562747.469425201,
       -5559747.066259384, -5556746.663093567,  -5553746.25992775,
       -5550745.856761932, -5547745.453596115, -5544745.050430298,
       -5541744.647264481,
       ...
        5538744.244098663,  5541744.647264481,  5544745.050430298,
        5547745.453596115,  5550745.856761932,   5553746.25992775,
        5556746.663093567,  5559747.066259384,  5562747.469425201,
        5565747.872591019],
      dtype='float64', name='x', length=3712))
  • orbital_parameters :
    {'projection_longitude': 0.0, 'projection_latitude': 0.0, 'projection_altitude': 35785831.0, 'satellite_nominal_longitude': 0.0, 'satellite_nominal_latitude': 0.0, 'satellite_actual_longitude': 0.14713798434861186, 'satellite_actual_latitude': -0.46057113067328453, 'satellite_actual_altitude': 35786185.500949726}
    time_parameters :
    {'nominal_start_time': datetime.datetime(2023, 1, 2, 12, 30), 'nominal_end_time': datetime.datetime(2023, 1, 2, 12, 45), 'observation_start_time': datetime.datetime(2023, 1, 2, 12, 30, 10, 117000), 'observation_end_time': datetime.datetime(2023, 1, 2, 12, 42, 43, 447000)}
    units :
    count
    wavelength :
    0.635 µm (0.56-0.71 µm)
    standard_name :
    counts
    platform_name :
    Meteosat-11
    sensor :
    seviri
    georef_offset_corrected :
    True
    start_time :
    2023-01-02 12:30:00
    end_time :
    2023-01-02 12:45:00
    reader :
    seviri_l1b_native
    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)
    name :
    VIS006
    resolution :
    3000.403165817
    calibration :
    counts
    modifiers :
    ()
    _satpy_id :
    DataID(name='VIS006', wavelength=WavelengthRange(min=0.56, central=0.635, max=0.71, unit='µm'), resolution=3000.403165817, calibration=<4>, modifiers=())
    ancillary_variables :
    []
In [10]:
print(scn['VIS006'].attrs['start_time'])
print(scn['VIS006'].attrs['end_time'])

2023-01-02 12:30:00
2023-01-02 12:45:00
In [11]:
print(scn)

<xarray.DataArray 'reshape-e4132bd22c06954dad3e0aad9b837671' (y: 3712, x: 3712)>
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(928, 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/18)
    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...
    time_parameters:          {'nominal_start_time': datetime.datetime(2023, ...
    units:                    count
    wavelength:               0.635 µm (0.56-0.71 µm)
    standard_name:            counts
    platform_name:            Meteosat-11
    ...                       ...
    name:                     VIS006
    resolution:               3000.403165817
    calibration:              counts
    modifiers:                ()
    _satpy_id:                DataID(name='VIS006', wavelength=WavelengthRang...
    ancillary_variables:      []
In [12]:
d = np.array(list(scn.values()))
In [13]:
np.shape(d)
Out[13]:
(1, 3712, 3712)
In [14]:
vis006 = scn['VIS006']
vis006.attrs['area']
Out[14]:
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)
In [15]:
vis006_meas = d[0]
vis006_meas_int = vis006_meas.astype(int)
print(vis006_meas_int[1600,1600])

115
In [16]:
print(scn.available_composite_names())

['airmass', 'ash', 'cloud_phase_distinction', 'cloud_phase_distinction_raw', '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_raw_with_night_ir', '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', 'rocket_plume_day', 'rocket_plume_night', 'snow', 'vis_sharpened_ir']
In [17]:
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)
In [18]:
rcNew.array2raster(vis006_meas.flatten())
rcNew.pix2value(ilwis.Pixel(1600,1600))
Out[18]:
115.0
In [19]:
plt.figure(figsize=(8, 8))
plt.imshow(vis006_meas, interpolation='nearest', cmap='Greys_r')
plt.show()
In [20]:
# 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
selection:  in 0.0366593 seconds
Out[20]:
288.0
In [21]:
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()
In [22]:
#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_1004799",GEOCS["_ANONYMOUS_1004799",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
In [23]:
b1s = ilwis.do('linearstretch',rcSelect, 1)
b1s = ilwis.do('setvaluerange', b1s, 0, 255, 1)

setvaluerange:  in 0.0859987 seconds
In [24]:
#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_1004799",GEOCS["_ANONYMOUS_1004799",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
In [25]:
b1s.pix2value(ilwis.Pixel(500,500))
Out[25]:
76.0
In [26]:
b1s.store('sub_stretched.mpr')
In [27]:
rcNew2 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew) # remove satpy nodata

mapcalc:  in 1.86548 seconds
In [28]:
rcNew2.store('vis006a.mpr')
In [29]:
scn.available_dataset_names()
Out[29]:
['HRV',
 'IR_016',
 'IR_039',
 'IR_087',
 'IR_097',
 'IR_108',
 'IR_120',
 'IR_134',
 'VIS006',
 'VIS008',
 'WV_062',
 'WV_073']
In [30]:
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])
In [31]:
from satpy.composites import GenericCompositor
comp = GenericCompositor('my_rgb')
my_rgb = comp((scn[r_band], scn[g_band], scn[b_band]))
In [32]:
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')
Out[32]:
<matplotlib.image.AxesImage at 0x7f57c01cdc50>
In [33]:
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.5361285
4.2887516
4.0803065
6.3270054
23.298351
20.061512
In [34]:
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)
In [35]:
rcNew2.array2raster(data)
In [36]:
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.327005386352539
23.298351287841797
20.061511993408203
In [37]:
rcNew2 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew2) # remove satpy nodata
rcNew2.store('msg_3_bands.mpl')

mapcalc:  in 7.04422 seconds
In [38]:
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_1004860",GEOCS["_ANONYMOUS_1004860",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
In [39]:
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)
In [40]:
b006.array2raster(vis006_rot.flatten())
print(b006.pix2value(ilwis.Pixel(1600,1600)))
print(b006.size())

6.327005386352539
Size(3712, 3712, 1)
In [41]:
b006 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', b006) # remove satpy nodata
b006.store('b006.mpr')

mapcalc:  in 1.97787 seconds
In [42]:
b006_sub = ilwis.do('select',b006, 'boundingbox(1300 250, 1649 499)')

selection:  in 0.00376804 seconds
In [43]:
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.372626781463623
Maximum: 42.50956344604492
Mean: 11.466866739063263
Median: 11.408899263637796
STD: 4.870460894488165
Count: 87500.0
In [44]:
print(b006_sub.size())

Size(350, 250, 1)
In [45]:
b006_sub.store('sub.mpr')
In [46]:
grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=20 40 40 20, gridsize=600 600, cornerofcorners=yes')
In [47]:
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
15.520933151245117
17.794700622558594
3.0646426677703857
0.0
resample:  in 0.298776 seconds
In [48]:
b006res.store('b006res.mpr')

Section on use of colour composites¶

In [50]:
scn4 = Scene(reader='seviri_l1b_native', filenames=fnames)

cir = 'colorized_ir_clouds'
scn4.load([cir], upper_right_corner='NE')
scn4.show(cir)
In [51]:
scn = Scene(reader='seviri_l1b_native', filenames=fnames)
scn.load(['colorized_ir_clouds'], upper_right_corner='NE')
scn['colorized_ir_clouds']
Out[51]:
<xarray.DataArray 'getitem-cf99d5b6e5785d5f5772ee9cbc5e8d08' (y: 3712, x: 3712)>
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(928, 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/21)
    orbital_parameters:       {'projection_longitude': 0.0, 'projection_latit...
    time_parameters:          {'nominal_start_time': datetime.datetime(2023, ...
    units:                    K
    wavelength:               10.8 µm (9.8-11.8 µm)
    standard_name:            colorized_ir_clouds
    platform_name:            Meteosat-11
    ...                       ...
    modifiers:                ()
    _satpy_id:                DataID(name='colorized_ir_clouds', resolution=3...
    ancillary_variables:      []
    optional_datasets:        []
    prerequisites:            [DataQuery(name='IR_108')]
    optional_prerequisites:   []
xarray.DataArray
'getitem-cf99d5b6e5785d5f5772ee9cbc5e8d08'
  • y: 3712
  • x: 3712
  • dask.array<chunksize=(928, 3712), meta=np.ndarray>
    Array Chunk
    Bytes 52.56 MiB 13.14 MiB
    Shape (3712, 3712) (928, 3712)
    Dask graph 4 chunks in 48 graph layers
    Data type float32 numpy.ndarray
    		3712 3712
    • acq_time
      (y)
      datetime64[ns]
      NaT NaT NaT NaT ... NaT NaT NaT NaT
      long_name :
      Mean scanline acquisition time
      array(['NaT', 'NaT', 'NaT', ..., 'NaT', 'NaT', 'NaT'],
      dtype='datetime64[ns]')
    • crs
      ()
      object
      PROJCRS["unknown",BASEGEOGCRS["u...
      array(<Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unk ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Geostationary Satellite (Sweep Y)
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich
, dtype=object)
    • y
      (y)
      float64
      5.569e+06 5.566e+06 ... -5.566e+06
      units :
      meter
      array([ 5568748.275757,  5565747.872591,  5562747.469425, ..., -5559747.066259,
       -5562747.469425, -5565747.872591])
    • x
      (x)
      float64
      -5.569e+06 -5.566e+06 ... 5.566e+06
      units :
      meter
      array([-5568748.275757, -5565747.872591, -5562747.469425, ...,  5559747.066259,
        5562747.469425,  5565747.872591])
    • y
      PandasIndex
      PandasIndex(Index([ 5568748.275756836,  5565747.872591019,  5562747.469425201,
        5559747.066259384,  5556746.663093567,   5553746.25992775,
        5550745.856761932,  5547745.453596115,  5544745.050430298,
        5541744.647264481,
       ...
       -5538744.244098663, -5541744.647264481, -5544745.050430298,
       -5547745.453596115, -5550745.856761932,  -5553746.25992775,
       -5556746.663093567, -5559747.066259384, -5562747.469425201,
       -5565747.872591019],
      dtype='float64', name='y', length=3712))
    • x
      PandasIndex
      PandasIndex(Index([-5568748.275756836, -5565747.872591019, -5562747.469425201,
       -5559747.066259384, -5556746.663093567,  -5553746.25992775,
       -5550745.856761932, -5547745.453596115, -5544745.050430298,
       -5541744.647264481,
       ...
        5538744.244098663,  5541744.647264481,  5544745.050430298,
        5547745.453596115,  5550745.856761932,   5553746.25992775,
        5556746.663093567,  5559747.066259384,  5562747.469425201,
        5565747.872591019],
      dtype='float64', name='x', length=3712))
  • orbital_parameters :
    {'projection_longitude': 0.0, 'projection_latitude': 0.0, 'projection_altitude': 35785831.0, 'satellite_nominal_longitude': 0.0, 'satellite_nominal_latitude': 0.0, 'satellite_actual_longitude': 0.14713798434861186, 'satellite_actual_latitude': -0.46057113067328453, 'satellite_actual_altitude': 35786185.500949726}
    time_parameters :
    {'nominal_start_time': datetime.datetime(2023, 1, 2, 12, 30), 'nominal_end_time': datetime.datetime(2023, 1, 2, 12, 45), 'observation_start_time': datetime.datetime(2023, 1, 2, 12, 30, 10, 117000), 'observation_end_time': datetime.datetime(2023, 1, 2, 12, 42, 43, 447000)}
    units :
    K
    wavelength :
    10.8 µm (9.8-11.8 µm)
    standard_name :
    colorized_ir_clouds
    platform_name :
    Meteosat-11
    sensor :
    seviri
    georef_offset_corrected :
    True
    start_time :
    2023-01-02 12:30:00
    end_time :
    2023-01-02 12:45:00
    reader :
    seviri_l1b_native
    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)
    name :
    colorized_ir_clouds
    resolution :
    3000.403165817
    calibration :
    brightness_temperature
    modifiers :
    ()
    _satpy_id :
    DataID(name='colorized_ir_clouds', resolution=3000.403165817)
    ancillary_variables :
    []
    optional_datasets :
    []
    prerequisites :
    [DataQuery(name='IR_108')]
    optional_prerequisites :
    []
In [52]:
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)
Out[52]:
<matplotlib.image.AxesImage at 0x7f56feec4490>
In [53]:
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()
In [54]:
data1 = np.array(data1_rot).flatten()
In [55]:
print(data1.ndim)
print(data1.shape)
print(data1.size)

1
(41336832,)
41336832
In [56]:
print(len(data1))

41336832
In [57]:
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)
In [58]:
rcNew3.array2raster(data1)
rcNew3 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew3) # remove satpy nodata

mapcalc:  in 6.32069 seconds
In [59]:
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
In [60]:
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_1004976",GEOCS["_ANONYMOUS_1004976",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
In [61]:
rcNew3.store('cir.mpl')

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

In [62]:
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']
In [63]:
#get listing of other functions
dir(crr)
Out[63]:
['__class__',
 '__contains__',
 '__delattr__',
 '__delitem__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__getitem__',
 '__getstate__',
 '__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_area_defs',
 '_compare_areas',
 '_compare_swath_defs',
 '_contained_sensor_names',
 '_copy_datasets_and_wishlist',
 '_create_reader_instances',
 '_datasets',
 '_dependency_tree',
 '_filter_loaded_datasets_from_trunk_nodes',
 '_gather_all_areas',
 '_generate_composite',
 '_generate_composites_from_loaded_datasets',
 '_generate_composites_nodes_from_loaded_datasets',
 '_get_finalized_destination_area',
 '_get_prereq_datasets',
 '_get_writer_by_ext',
 '_ipython_key_completions_',
 '_load_datasets_by_readers',
 '_prepare_resampler',
 '_read_dataset_nodes_from_storage',
 '_read_datasets_from_storage',
 '_reader_times',
 '_readers',
 '_reduce_data',
 '_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',
 'chunk',
 'coarsest_area',
 'compute',
 'copy',
 'crop',
 'end_time',
 'finest_area',
 'generate_possible_composites',
 'get',
 'images',
 'iter_by_area',
 'keys',
 'load',
 'max_area',
 'min_area',
 'missing_datasets',
 'persist',
 'resample',
 'save_dataset',
 'save_datasets',
 'sensor_names',
 'show',
 'slice',
 'start_time',
 'to_geoviews',
 'to_xarray',
 'to_xarray_dataset',
 'unload',
 'values',
 'wishlist']
In [64]:
crr.load(['convective_precipitation_hourly_accumulation', 'convective_rain_rate'])
In [65]:
#crr.load(['convective_precipitation_hourly_accumulation'])
In [66]:
print(crr)

<xarray.DataArray 'add-5cae7ab9b61655233650b20852652db0' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=uint8, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/24)
    start_time:              2023-01-17 05:15:00
    end_time:                2023-01-17 05:27:22
    orbital_parameters:      {'satellite_nominal_altitude': 35785863.0, 'sate...
    long_name:               NWC GEO CRR Convective Rainfall Rate Class
    valid_range:             [ 0 11]
    _FillValue:              255
    ...                      ...
    modifiers:               ()
    _satpy_id:               DataID(name='convective_rain_rate', resolution=3...
    optional_datasets:       []
    standard_name:           convective_rain_rate
    prerequisites:           ['crr']
    optional_prerequisites:  []
<xarray.DataArray 'add-5bca5bb75ee1e01c90ada14bce8d1318' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/22)
    start_time:              2023-01-17 05:15:00
    end_time:                2023-01-17 05:27:22
    orbital_parameters:      {'satellite_nominal_altitude': 35785863.0, 'sate...
    long_name:               NWC GEO CRR Convective Hourly Rainfall Accumulation
    units:                   mm
    valid_range:             [ 0. 50.]
    ...                      ...
    modifiers:               ()
    _satpy_id:               DataID(name='convective_precipitation_hourly_acc...
    optional_datasets:       []
    standard_name:           convective_precipitation_hourly_accumulation
    prerequisites:           ['crr_accum']
    optional_prerequisites:  []
In [67]:
crr.show('convective_rain_rate')
Out[67]:
In [68]:
crr_load = crr['convective_rain_rate']
In [69]:
print(crr_load)
print(crr_load.shape)

<xarray.DataArray 'add-5cae7ab9b61655233650b20852652db0' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=uint8, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/24)
    start_time:              2023-01-17 05:15:00
    end_time:                2023-01-17 05:27:22
    orbital_parameters:      {'satellite_nominal_altitude': 35785863.0, 'sate...
    long_name:               NWC GEO CRR Convective Rainfall Rate Class
    valid_range:             [ 0 11]
    _FillValue:              255
    ...                      ...
    modifiers:               ()
    _satpy_id:               DataID(name='convective_rain_rate', resolution=3...
    optional_datasets:       []
    standard_name:           convective_rain_rate
    prerequisites:           ['crr']
    optional_prerequisites:  []
(3712, 3712)
In [70]:
crr_meas = crr_load.values
In [71]:
print(crr_meas)
print(crr_meas.shape)

[[255 255 255 ... 255 255 255]
 [255 255 255 ... 255 255 255]
 [255 255 255 ... 255 255 255]
 ...
 [255 255 255 ... 255 255 255]
 [255 255 255 ... 255 255 255]
 [255 255 255 ... 255 255 255]]
(3712, 3712)
In [72]:
plt.imshow(crr_meas, interpolation='nearest', vmin=0, vmax=11, cmap='jet')
plt.show()
In [73]:
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,1))
rcNew.setGeoReference(grf)
rcNew.setDataDef(dfNum)
In [74]:
rcNew.array2raster(crr_meas.flatten())
rcNew = ilwis.do('mapcalc', 'iff(@1<255,@1,?)', rcNew) # remove background - assign to no data

mapcalc:  in 1.69602 seconds
In [75]:
print(rcNew.coord2value(ilwis.Coordinate(4896095.17,-1425952.85))) 

11.0
In [76]:
rcNew.store('crr_class.mpr')
In [77]:
crr.load(['crr_accum'])
In [78]:
print(crr)

<xarray.DataArray 'add-5cae7ab9b61655233650b20852652db0' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=uint8, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/24)
    start_time:              2023-01-17 05:15:00
    end_time:                2023-01-17 05:27:22
    orbital_parameters:      {'satellite_nominal_altitude': 35785863.0, 'sate...
    long_name:               NWC GEO CRR Convective Rainfall Rate Class
    valid_range:             [ 0 11]
    _FillValue:              255
    ...                      ...
    modifiers:               ()
    _satpy_id:               DataID(name='convective_rain_rate', resolution=3...
    optional_datasets:       []
    standard_name:           convective_rain_rate
    prerequisites:           ['crr']
    optional_prerequisites:  []
<xarray.DataArray 'add-5bca5bb75ee1e01c90ada14bce8d1318' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/22)
    start_time:              2023-01-17 05:15:00
    end_time:                2023-01-17 05:27:22
    orbital_parameters:      {'satellite_nominal_altitude': 35785863.0, 'sate...
    long_name:               NWC GEO CRR Convective Hourly Rainfall Accumulation
    units:                   mm
    valid_range:             [ 0. 50.]
    ...                      ...
    modifiers:               ()
    _satpy_id:               DataID(name='convective_precipitation_hourly_acc...
    optional_datasets:       []
    standard_name:           convective_precipitation_hourly_accumulation
    prerequisites:           ['crr_accum']
    optional_prerequisites:  []
<xarray.DataArray 'crr_accum' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=float32, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/18)
    start_time:           2023-01-17 05:15:00
    end_time:             2023-01-17 05:27:22
    orbital_parameters:   {'satellite_nominal_altitude': 35785863.0, 'satelli...
    long_name:            NWC GEO CRR Convective Hourly Rainfall Accumulation
    units:                mm
    valid_range:          [ 0. 50.]
    ...                   ...
    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...
In [79]:
crr.show('crr_accum')
Out[79]:
In [80]:
crr_meas = crr['crr_accum'].values
In [81]:
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)
In [82]:
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))) 

mapcalc:  in 1.92344 seconds22.399999618530273
In [83]:
rcNew.store('crr_value.mpr')
In [84]:
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']
In [85]:
rdt.load(['MapCellCatType', 'MapCellCatType_pal'])
In [86]:
rdt.load(['rdt_cell_type'])
In [87]:
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) float32 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/21)
    start_time:           2023-01-17 03:45:00
    end_time:             2023-01-17 03:57:22
    orbital_parameters:   {'satellite_nominal_altitude': 35785863.0, 'satelli...
    long_name:            NWC GEO RDT-CW Type and phase of significant cells
    valid_range:          [ 0 21]
    _FillValue:           255
    ...                   ...
    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/16)
    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 'add-d0abc54d16ee05d9bf3ef6ab05f3ca7d' (y: 3712, x: 3712)>
dask.array<add, shape=(3712, 3712), dtype=uint8, chunksize=(3712, 3712), chunktype=numpy.ndarray>
Coordinates:
  * y        (y) float32 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/25)
    start_time:              2023-01-17 03:45:00
    end_time:                2023-01-17 03:57:22
    orbital_parameters:      {'satellite_nominal_altitude': 35785863.0, 'sate...
    long_name:               NWC GEO RDT-CW Type and phase of significant cells
    valid_range:             [ 0 21]
    _FillValue:              255
    ...                      ...
    modifiers:               ()
    _satpy_id:               DataID(name='rdt_cell_type', resolution=3000)
    optional_datasets:       []
    standard_name:           rdt_cell_type
    prerequisites:           ['MapCellCatType']
    optional_prerequisites:  []
In [88]:
rdt.show('rdt_cell_type')
Out[88]:
In [89]:
rdt_meas = rdt['rdt_cell_type'].values
In [90]:
print(rdt_meas.shape)

(3712, 3712)
In [91]:
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)
In [92]:
rcNew.array2raster(rdt_meas.flatten())
rcNew = ilwis.do('mapcalc', 'iff(@1>0,@1,?)', rcNew) # remove no rainall (value = 0) and assign it to nodata

mapcalc:  in 1.89936 seconds
In [93]:
rcNew.store('rdt_class.mpr')

Radipdly Developing Thunderstorms classes description:

  • 0=Non_convective
  • 1=Convective_triggering
  • 2=Convective_triggering_from_split
  • 3=Convective_growing
  • 4=Convective_mature
  • 5=vershootingTop_mature
  • 6=Convective_decaying
  • 7=Electric_triggering
  • 8=Electric_triggering_from_split
  • 9=Electric_growing
  • 10=Electric_mature
  • 11=Electric_decaying
  • 12=HighRainRate_triggering
  • 13=HighRainRate_triggering_from_split
  • 14=HighRainRate_growing
  • 15=HighRainRate_mature
  • 16=HighRainRate_decaying
  • 17=HighSeverity_triggering
  • 18=HighSeverity_triggering_from_split
  • 19=HighSeverity_growing
  • 20=HighSeverity_mature
  • 21=HighSeverity_decaying

Review all files created using ILWIS386