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]:
import ilwis
import datetime
import shutil
from IPython.core.display import HTML
import os
import os.path
import sys
import glob
import numpy as np
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')
In [2]:
#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)
/home/ilwispy/notebooks/ilwispy_tutorial/Intro_Satpy2ILWISPy
In [3]:
ilwis.version()
Out[3]:
'1.0 build 20250227'
In [4]:
# read unzipped file, here extension *.nat
fnames = glob.glob(work_dir+'/*.nat')
print(fnames)
scn = Scene(reader='seviri_l1b_native', filenames=fnames)
['/home/ilwispy/notebooks/ilwispy_tutorial/Intro_Satpy2ILWISPy/MSG4-SEVI-MSG15-0100-NA-20230102124243.447000000Z-NA.nat']
In [5]:
scn.available_dataset_names()
Out[5]:
['HRV', 'IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073']
In [6]:
#get list of functions available
dir(scn)
Out[6]:
['__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_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_hvplot', 'to_xarray', 'to_xarray_dataset', 'unload', 'values', 'wishlist']
In [7]:
#note calibration = counts
scn.load(['VIS006'], upper_right_corner='NE', calibration='counts')
scn.keys()
Out[7]:
[DataID(name='VIS006', wavelength=WavelengthRange(min=0.56, central=0.635, max=0.71, unit='µm'), resolution=3000.403165817, calibration=<calibration.counts>, modifiers=())]
In [8]:
scn['VIS006']
Out[8]:
<xarray.DataArray 'reshape-5b90d6381d47853c787158fa86167ddc' (y: 3712, x: 3712)> Size: 55MB
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(460, 3712), chunktype=numpy.ndarray>
Coordinates:
acq_time (y) datetime64[ns] 30kB NaT NaT NaT NaT NaT ... NaT NaT NaT NaT
crs object 8B PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknow...
* y (y) float64 30kB 5.569e+06 5.566e+06 ... -5.563e+06 -5.566e+06
* x (x) float64 30kB -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
Attributes: (12/18)
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: []In [9]:
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 [10]:
print(scn)
<xarray.DataArray 'reshape-5b90d6381d47853c787158fa86167ddc' (y: 3712, x: 3712)> Size: 55MB
dask.array<getitem, shape=(3712, 3712), dtype=float32, chunksize=(460, 3712), chunktype=numpy.ndarray>
Coordinates:
acq_time (y) datetime64[ns] 30kB NaT NaT NaT NaT NaT ... NaT NaT NaT NaT
crs object 8B PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknow...
* y (y) float64 30kB 5.569e+06 5.566e+06 ... -5.563e+06 -5.566e+06
* x (x) float64 30kB -5.569e+06 -5.566e+06 ... 5.563e+06 5.566e+06
Attributes: (12/18)
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: []
In [11]:
d = np.array(list(scn.values()))
In [12]:
np.shape(d)
Out[12]:
(1, 3712, 3712)
In [13]:
vis006 = scn['VIS006']
vis006.attrs['area']
Out[13]:
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 [14]:
vis006_meas = d[0]
vis006_meas_int = vis006_meas.astype(int)
print(vis006_meas_int[1600,1600])
115
In [15]:
print(scn.available_composite_names())
['24h_microphysics', '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 [16]:
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 [17]:
rcNew.array2raster(vis006_meas.flatten())
rcNew.pix2value(ilwis.Pixel(1600,1600))
Out[17]:
115.0
In [18]:
plt.figure(figsize=(8, 8))
plt.imshow(vis006_meas, interpolation='nearest', cmap='Greys_r')
plt.show()
In [19]:
# 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
Out[19]:
288.0
In [20]:
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 [21]:
#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_1004730",GEOCS["_ANONYMOUS_1004730",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 [22]:
b1s = ilwis.do('linearstretch',rcSelect, 1)
b1s = ilwis.do('setvaluerange', b1s, 0, 255, 1)
In [23]:
#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_1004730",GEOCS["_ANONYMOUS_1004730",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 [24]:
b1s.pix2value(ilwis.Pixel(500,500))
Out[24]:
76.0
In [25]:
b1s.store('sub_stretched.mpr')
In [26]:
rcNew2 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew) # remove satpy nodata
In [27]:
rcNew2.store('vis006a.mpr')
In [28]:
scn.available_dataset_names()
Out[28]:
['HRV', 'IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073']
In [29]:
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 [30]:
from satpy.composites import GenericCompositor
comp = GenericCompositor('my_rgb')
my_rgb = comp((scn[r_band], scn[g_band], scn[b_band]))
In [31]:
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[31]:
<matplotlib.image.AxesImage at 0x7cbc7c44b130>
In [32]:
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.536129 4.2887516 4.080307 6.327006 23.298353 20.061514
In [33]:
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 [34]:
rcNew2.array2raster(data)
In [35]:
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.327005863189697 23.29835319519043 20.061513900756836
In [36]:
rcNew2 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', rcNew2) # remove satpy nodata
rcNew2.store('msg_3_bands.mpl')
In [37]:
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_1004809",GEOCS["_ANONYMOUS_1004809",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 [38]:
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 [39]:
b006.array2raster(vis006_rot.flatten())
print(b006.pix2value(ilwis.Pixel(1600,1600)))
print(b006.size())
6.327005863189697 Size(3712, 3712, 1)
In [40]:
b006 = ilwis.do('mapcalc', 'iff(@1>=0,@1,?)', b006) # remove satpy nodata
b006.store('b006.mpr')
In [41]:
b006_sub = ilwis.do('select',b006, 'boundingbox(1300 250, 1649 499)')
In [42]:
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.372627019882202 Maximum: 42.50956344604492 Mean: 11.466867353354862 Median: 11.40890032358629 STD: 4.870461139023392 Count: 87500.0
In [43]:
print(b006_sub.size())
Size(350, 250, 1)
In [44]:
b006_sub.store('sub.mpr')
In [45]:
grf_LL= ilwis.GeoReference('code=georef:type=corners, csy=epsg:4326, envelope=20 40 40 20, gridsize=600 600, cornerofcorners=yes')
In [46]:
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.520934104919434 17.794702529907227 3.064642906188965 0.0
In [47]:
b006res.store('b006res.mpr')
Section on use of colour composites¶
In [48]:
scn4 = Scene(reader='seviri_l1b_native', filenames=fnames)
cir = 'colorized_ir_clouds'
scn4.load([cir], upper_right_corner='NE')
scn4.show(cir)
Out[48]:
