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: []
|
array(['NaT', 'NaT', 'NaT', ..., 'NaT', 'NaT', 'NaT'], dtype='datetime64[ns]')
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)
array([ 5568748.275757, 5565747.872591, 5562747.469425, ..., -5559747.066259, -5562747.469425, -5565747.872591])
array([-5568748.275757, -5565747.872591, -5562747.469425, ..., 5559747.066259, 5562747.469425, 5565747.872591])
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')
scn4 = Scene(reader='seviri_l1b_native', filenames=fnames)
cir = 'colorized_ir_clouds'
scn4.load([cir], upper_right_corner='NE')
scn4.show(cir)