From 6bf44acfcbfc6d4ee2c67c57cc03fdb147eb73ad Mon Sep 17 00:00:00 2001 From: ArchiPEL v3 Date: Thu, 26 Sep 2024 14:37:30 +0000 Subject: [PATCH 1/9] A new reader for the netcdf sent to Icare by Meteo-France. --- satpy/etc/readers/goesr_netcdficare.yaml | 309 +++++++++ satpy/etc/readers/hima_netcdficare.yaml | 261 +++++++ satpy/etc/readers/msg_netcdficare.yaml | 174 +++++ satpy/etc/readers/mtg_netcdficare.yaml | 335 +++++++++ satpy/readers/geos_netcdficare.py | 826 +++++++++++++++++++++++ 5 files changed, 1905 insertions(+) create mode 100644 satpy/etc/readers/goesr_netcdficare.yaml create mode 100644 satpy/etc/readers/hima_netcdficare.yaml create mode 100644 satpy/etc/readers/msg_netcdficare.yaml create mode 100644 satpy/etc/readers/mtg_netcdficare.yaml create mode 100644 satpy/readers/geos_netcdficare.py diff --git a/satpy/etc/readers/goesr_netcdficare.yaml b/satpy/etc/readers/goesr_netcdficare.yaml new file mode 100644 index 0000000000..7118e0dd07 --- /dev/null +++ b/satpy/etc/readers/goesr_netcdficare.yaml @@ -0,0 +1,309 @@ +# References: +# - GOES Level 1.5 Image Data Format Description +# - Radiometric Calibration of GOESR ABI Level 1.5 Image Data in Equivalent +# Spectral Blackbody Radiance +# Netcdf built by Icare and Meteo France. Stored at Icare. + +reader: + name: goesr_netcdficare + short_name: goesr_NETCDF_ICARE + long_name: METEOFRANCE GEOSTATIONARY NETCDF BUILD for ICARE (Lille) + description: A reader for L1b NETCDF for all GEOSTATIONNARY retrieved from the ICARE service. + status: Defunct + supports_fsspec: false + sensors: [abi] + default_channels: [VIS_004, VIS_006, VIS_008, VIS_014, VIS_016, VIS_022, IR_039, IR_062, IR_069, IR_073, IR_085, IR_096, IR_103, IR_112, IR_123, IR_133] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + + GOES500m : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Emultic500mNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Wmultic500mNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Emultic500mNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc', + 'Wmultic500mNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc'] + + GOES1km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Emultic1kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Wmultic1kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Emultic1kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc', + 'Wmultic1kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc'] + + GOES2km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Emultic2kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Wmultic2kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Emultic2kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc', + 'Wmultic2kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc'] + + +datasets: + + VIS_004 : + name: C01 + wavelength: [0.450, 0.470, 0.490] + resolution: + 1000: { file_type: GOES1km } + 2000: { file_type: GOES2km } + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: "1" + file_type: [GOES2km, GOES1km] + + VIS_006 : + name: C02 + wavelength: [0.590, 0.640, 0.690] + resolution: + 500: { file_type: GOES500m } + 2000: { file_type: GOES2km } + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: "1" + file_type: [GOES2km, GOES500m] + + VIS_008 : + name: C03 + wavelength: [0.8455, 0.865, 0.8845] + resolution: + 1000: { file_type: GOES1km } + 2000: { file_type: GOES2km } + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: "1" + file_type: [GOES2km, GOES1km] + + VIS_014 : + name: C04 + wavelength: [1.3705, 1.378, 1.3855] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + VIS_016 : + name: C05 + wavelength: [1.580, 1.610, 1.640] + resolution: + 1000: { file_type: GOES1km } + 2000: { file_type: GOES2km } + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: "1" + file_type: [GOES2km, GOES1km] + + VIS_022 : + name: C06 + wavelength: [2.225, 2.250, 2.275] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavelength + units: W m-2 um-1 sr-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_039 : + name: C07 + wavelength: [3.80, 3.90, 4.00] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_062 : + name: C08 + wavelength: [5.770, 6.185, 6.600] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_069 : + name: C09 + wavelength: [6.75, 6.95, 7.15] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_073 : + name: C10 + wavelength: [7.24, 7.34, 7.44] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_085 : + name: C11 + wavelength: [8.30, 8.50, 8.70] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_096 : + name: C12 + wavelength: [9.42, 9.61, 9.80] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_103 : + name: C13 + wavelength: [10.10, 10.35, 10.60] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_112 : + name: C14 + wavelength: [10.80, 11.20, 11.60] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_123 : + name: C15 + wavelength: [11.80, 12.30, 12.80] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + IR_133 : + name: C16 + wavelength: [13.00, 13.30, 13.60] + resolution: 2000 + calibration: + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: K + counts: + standard_name: counts + units: "1" + file_type: GOES2km + + + diff --git a/satpy/etc/readers/hima_netcdficare.yaml b/satpy/etc/readers/hima_netcdficare.yaml new file mode 100644 index 0000000000..cdd3bcfd96 --- /dev/null +++ b/satpy/etc/readers/hima_netcdficare.yaml @@ -0,0 +1,261 @@ +# References: +# - Himawari Level 1.5 Image Data Format Description +# - Radiometric Calibration of Himawari Level 1.5 Image Data in Equivalent +# Spectral Blackbody Radiance +# Netcdf built by Icare and Meteo France. Stored at Icare. + +reader: + name: hima_netcdficare + short_name: hima_NETCDF_ICARE + long_name: METEOFRANCE GEOSTATIONARY NETCDF BUILD for ICARE (Lille) + description: A reader for L1b NETCDF for all GEOSTATIONNARY retrieved from the ICARE service. + status: Defunct + supports_fsspec: false + sensors: [ahi] + default_channels: [VIS004, VIS005, VIS006, VIS008, IR_016, IR_022, IR_038, WV_062, WV_069, WV_073, IR_085, IR_096, IR_104, IR_112, IR_123, IR_132] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + + HIMA500m : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Jmultic500mNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Jmultic500mNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc'] + + HIMA1km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Jmultic1kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Jmultic1kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc'] + + HIMA2km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Jmultic2kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d%H%M}.nc', + 'Jmultic2kmNC4_{platform_shortname:6s}_{start_time:%Y%m%d_%H%M}.nc'] + +datasets: + + VIS004: + name: B01 + sensor: ahi + wavelength: [0.45,0.47,0.49] + resolution: [2000] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + VIS005: + name: B02 + sensor: ahi + wavelength: [0.49,0.51,0.53] + resolution: [2000] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + VIS006: + name: B03 + sensor: ahi + wavelength: [0.62,0.64,0.66] + resolution: [500, 2000] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: 1 + file_type: [HIMA500m, HIMA2km] + + VIS008 : + name: B04 + sensor: ahi + wavelength: [0.83, 0.85, 0.87] + resolution: [1000, 2000] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: 1 + file_type: [HIMA1km, HIMA2km] + + IR_016 : + name: B05 + sensor: ahi + wavelength: [1.5, 1.6, 1.7] + resolution: 2000 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_022 : + name: B06 + sensor: ahi + wavelength: [2.2, 2.3, 2.4] + resolution: 2000 + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_038 : + name: B07 + resolution: 2000 + sensor: ahi + wavelength: [3.7, 3.9, 4.1] + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + # FUTURE: Split this in to multiple resolutions so each can be loaded + file_type: [HIMA2km] + + WV_062 : + name: B08 + sensor: ahi + wavelength: [6.0, 6.2, 6.4] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + WV_069 : + name: B09 + sensor: ahi + wavelength: [6.7, 6.9, 7.1] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + WV_073 : + name: B10 + sensor: ahi + wavelength: [7.1, 7.3, 7.5] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_085 : + name: B11 + sensor: ahi + wavelength: [8.4, 8.6, 8.8] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_096 : + name: B12 + sensor: ahi + wavelength: [9.4, 9.6, 9.8] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_104 : + name: B13 + sensor: ahi + wavelength: [10.2, 10.4, 10.6] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_112 : + name: B14 + sensor: ahi + wavelength: [11.0, 11.2, 11.4] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_123 : + name: B15 + sensor: ahi + wavelength: [12.2, 12.4, 12.6] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + + IR_132 : + name: B16 + sensor: ahi + wavelength: [13.1, 13.3, 13.5] + resolution: 2000 + calibration: + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + counts: + standard_name: counts + units: 1 + file_type: [HIMA2km] + diff --git a/satpy/etc/readers/msg_netcdficare.yaml b/satpy/etc/readers/msg_netcdficare.yaml new file mode 100644 index 0000000000..01fa4245cb --- /dev/null +++ b/satpy/etc/readers/msg_netcdficare.yaml @@ -0,0 +1,174 @@ +# References: +# - MSG Level 1.5 Image Data Format Description +# - Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent +# Spectral Blackbody Radiance + +reader: + name: msg_netcdficare + short_name: MSG L1B NETCDF ICARE METEO-FRANCE + long_name: METEOFRANCE GEOSTATIONARY NETCDF BUILD for ICARE (Lille) + description: A reader for L1b NETCDF for all GEOSTATIONNARY retrieved from the ICARE service. + status: Defunct + supports_fsspec: false + sensors: [seviri] + default_channels: [HRV, IR_016, IR_039, IR_087, IR_097, IR_108, IR_120, IR_134, VIS006, VIS008, WV_062, WV_073] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + + MSG3km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Mmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Imultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mrsmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Imultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Mrsmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc'] + + # platform_shortname:5s : msg01, ..., msg04 + + MSG1km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: [ 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc'] + # platform_shortname:5s : msg01, ..., msg04 + +datasets: + + HRV: + name: HRV + sensor: seviri + resolution: 1000.134348869 + wavelength: [0.5, 0.7, 0.9] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: MSG1km + + IR_016: + name: IR_016 + sensor: seviri + resolution: 3000.403165817 + wavelength: [1.5, 1.64, 1.78] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: MSG3km + + IR_039: + name: IR_039 + sensor: seviri + resolution: 3000.403165817 + wavelength: [3.48, 3.92, 4.36] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + IR_087: + name: IR_087 + sensor: seviri + resolution: 3000.403165817 + wavelength: [8.3, 8.7, 9.1] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + IR_097: + name: IR_097 + sensor: seviri + resolution: 3000.403165817 + wavelength: [9.38, 9.66, 9.94] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + IR_108: + name: IR_108 + sensor: seviri + resolution: 3000.403165817 + wavelength: [9.8, 10.8, 11.8] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + IR_120: + name: IR_120 + sensor: seviri + resolution: 3000.403165817 + wavelength: [11.0, 12.0, 13.0] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + IR_134: + name: IR_134 + sensor: seviri + resolution: 3000.403165817 + wavelength: [12.4, 13.4, 14.4] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + VIS006: + name: VIS006 + sensor: seviri + resolution: 3000.403165817 + wavelength: [0.56, 0.635, 0.71] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: MSG3km + + VIS008: + name: VIS008 + sensor: seviri + resolution: 3000.403165817 + wavelength: [0.74, 0.81, 0.88] + calibration: + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: MSG3km + + WV_062: + name: WV_062 + sensor: seviri + resolution: 3000.403165817 + wavelength: [5.35, 6.25, 7.15] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + WV_073: + name: WV_073 + sensor: seviri + resolution: 3000.403165817 + wavelength: [6.85, 7.35, 7.85] + calibration: + brightness_temperature: + standard_name: brightness_temperature + units: "K" + file_type: MSG3km + + diff --git a/satpy/etc/readers/mtg_netcdficare.yaml b/satpy/etc/readers/mtg_netcdficare.yaml new file mode 100644 index 0000000000..2562feafd7 --- /dev/null +++ b/satpy/etc/readers/mtg_netcdficare.yaml @@ -0,0 +1,335 @@ +# References: +# - MTG Level 1.5 Image Data Format Description +# - Radiometric Calibration of MSG SEVIRI Level 1.5 Image Data in Equivalent +# Spectral Blackbody Radiance +# Netcdf built by Icare and Meteo France. Stored at Icare. + +reader: + name: mtg_netcdficare + short_name: mtg_NETCDF_ICARE + long_name: METEOFRANCE GEOSTATIONARY NETCDF BUILD for ICARE (Lille) + description: A reader for L1b NETCDF for all GEOSTATIONNARY retrieved from the ICARE service. + status: Defunct + supports_fsspec: false + sensors: [fci] + default_channels: [VIS004, VIS005, VIS006, VIS008, VIS009, IR_013, IR_016, IR_022, IR_038, WV_063, WV_073, IR_087, IR_097, IR_105, IR_123, IR_133] + reader: !!python/name:satpy.readers.yaml_reader.FileYAMLReader + +file_types: + + MTG500m : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Multic500m_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Multic500m_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M%S}.nc', + 'Mmultic500mNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mmultic500mNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M%S}.nc'] + + MTG1km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Multic1km_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Multic1km_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M%S}.nc', + 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M%S}.nc'] + + MTG2km : + file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE + file_patterns: ['Multic2km_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Multic2km_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M%S}.nc', + 'Mmultic2kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mmultic2kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M%S}.nc'] + +datasets: + + VIS004: + name: vis_04 + sensor: fci + wavelength: [ 0.384, 0.444, 0.504 ] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: [MTG2km, MTG1km] + # file_type: abi_ahi_seviri_fci_netcdficare + # file_type: [mtg2km__netcdficare, mtg1km__netcdficare] + + VIS005: + name: vis_05 + sensor: fci + wavelength: [0.470, 0.510, 0.550] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: [MTG2km, MTG1km] + + VIS006: + name: vis_06 + sensor: fci + wavelength: [0.590, 0.640, 0.690] + resolution: + 500: { file_type: MTG500m } + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: [MTG2km, MTG1km, MTG500m] + + VIS008 : + name: vis_08 + sensor: fci + wavelength: [0.815, 0.865, 0.915] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: [MTG2km, MTG1km] + + VIS009 : + name: vis_09 + sensor: fci + wavelength: [0.894, 0.914, 0.934] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + + IR_013 : + name: nir_13 + sensor: fci + wavelength: [1.350, 1.380, 1.410] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + + IR_016 : + name: nir_16 + sensor: fci + wavelength: [1.560, 1.610, 1.660] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + file_type: [MTG2km, MTG1km] + + IR_022 : + name: nir_22 + sensor: fci + wavelength: [2.200, 2.250, 2.300] + resolution: + 500: { file_type: MTG500m } + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + reflectance: + standard_name: toa_bidirectional_reflectance + units: "%" + + IR_038 : + name: ir_38 + sensor: fci + wavelength: [3.400, 3.800, 4.200] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + file_type: [MTG2km, MTG1km] + + WV_063 : + name: wv_63 + sensor: fci + wavelength: [5.300, 6.300, 7.300] + resolution: + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + WV_073 : + name: wv_73 + sensor: fci + wavelength: [6.850, 7.350, 7.850] + resolution: + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + IR_087 : + name: ir_87 + sensor: fci + wavelength: [8.300, 8.700, 9.100] + resolution: + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + IR_097 : + name: ir_97 + sensor: fci + wavelength: [9.360, 9.660, 9.960] + resolution: + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + IR_105 : + name: ir_105 + sensor: fci + wavelength: [9.800, 10.500, 11.200] + resolution: + 1000: { file_type: MTG1km } + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + IR_123 : + name: ir_123 + sensor: fci + wavelength: [11.800, 12.300, 12.800] + resolution: + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + IR_133 : + name: ir_133 + sensor: fci + wavelength: [12.700, 13.300, 13.900] + resolution: + 2000: { file_type: MTG2km } + calibration: + counts: + standard_name: counts + units: "count" + radiance: + standard_name: toa_outgoing_radiance_per_unit_wavenumber + units: mW m-2 sr-1 (cm-1)-1 + brightness_temperature: + standard_name: toa_brightness_temperature + units: "K" + + diff --git a/satpy/readers/geos_netcdficare.py b/satpy/readers/geos_netcdficare.py new file mode 100644 index 0000000000..31bae413b9 --- /dev/null +++ b/satpy/readers/geos_netcdficare.py @@ -0,0 +1,826 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2019 Satpy developers +# +# This file is part of satpy. Written by Meteo France in august 2024. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with +# satpy. If not, see . + +"""Interface to GEOSTATIONNARY L1B NETCDF data from ICARE (Lille). + +Introduction +------------ + +The ``geos_netcdficare`` reader reads some geostationnary netcdf build by Meteo France and stored at Icare. +The brightness tempeture and albedo are calibrated. + +That has been stored by the ICARE Data and Services Center +Data can be accessed via: http://www.icare.univ-lille1.fr + +This reader concerns the following netcdfs : + +. msg with a longitude near 0° : +Mmultic3kmNC4_msg03_20231113_111500.nc Mmultic1kmNC4_msg03_20231113_111500.nc + +. Msg rapid scan with a longitude near 9.5° : +Mrsmultic3kmNC4_msg03_20231113_111500.nc Mrsmultic1kmNC4_msg03_20231113_111500.nc + +. Msg with a longitude near 42° : +Imultic3kmNC4_msg03_20231113_111500.nc Imultic1kmNC4_msg03_20231113_111500.nc + +. Himawari : +Jmultic2kmNC4_hima09_20231113_111500.nc Jmultic1kmNC4_hima09_20231113_111500.nc Jmultic500mNC4_hima09_20231113_111500.nc + +. Goesr near -137° : +Wmultic2kmNC4_goes16_202406281000.nc. The better resolution are not built at Lannion, only at Tahiti. + +. Goesr in -75° : +Emultic2kmNC4_goes16_202406281000.nc Emultic1kmNC4_goes16_202406281000.nc Emultic500mNC4_goes16_202406281000.nc + +. Mtg : +Mmultic2km_mtgi1_20240104_090000.nc Mmultic1km_mtgi1_20240104_090000.nc Mmultic500m_mtgi1_20240104_090000.nc + + +Example: +-------- +Here is an example how to read the data in satpy: + +.. code-block:: python + + from satpy import Scene + import glob + + filenames = glob.glob('data/*2019-03-01T12-00-00*.hdf') + scn = Scene(filenames = filenames, reader = 'hima_netcdficare') + scn.load(['true_color']) # scn.load(['VIS006']) + + my_area = AreaDefinition('my_area', 'zone', 'my_area', + '+proj=latlong +lon_0=0 +a=6378169 +b=6356583 +h=35785831 +x_0=0 +y_0=0 +pm=0', + 8500, 4000, + [-180., -80., 180., 80], + nprocs=16) + + natscn = scn.resample(my_area, resampler='nearest') + natscn.save_dataset(composite_name, filename = filename_image_out) + +Output: + +.. code-block:: none + + + dask.array + Coordinates: + crs object +proj=geos +a=6378169.0 +b=6356583.8 +lon_0=0.0 +h=35785831.0 +units=m +type=crs + * y (y) float64 5.566e+06 5.563e+06 5.56e+06 ... -5.566e+06 -5.569e+06 + * x (x) float64 -5.566e+06 -5.563e+06 -5.56e+06 ... 5.566e+06 5.569e+06 + Attributes: + start_time: 2004-12-29 12:15:00 + end_time: 2004-12-29 12:27:44 + area: Area ID: geosmsg\nDescription: MSG/SEVIRI low resol... + name: IR_108 + resolution: 3000.403165817 + calibration: brightness_temperature + polarization: None + level: None + modifiers: () + ancillary_variables: [] + + +""" + +import datetime as dt + +import numpy as np + +from satpy.readers._geos_area import get_area_definition, get_area_extent + +from netCDF4 import Dataset +from xarray import DataArray +import xarray as xr + +from satpy.readers.file_handlers import BaseFileHandler + +from satpy._compat import cached_property +# from satpy.utils import get_dask_chunk_size_in_bytes +import dask +from satpy.readers import open_file_or_filename +import math + +class NETCDF_ICARE(BaseFileHandler) : + # Cf readers/file_handlers.py. + + def __init__(self, filename, filename_info, filetype_info) : + """Init the file handler.""" + + super().__init__(filename, filename_info, filetype_info) # cache_var_size=0, cache_handle = False) + # self.ncfileIn = Dataset(filename, mode = "r", clobber=True, diskless=False, persist=False, format='NETCDF4') + + # self.filename = filename + # self.filetype_info = filetype_info + + # chunk_bytes = self._chunk_bytes_for_resolution() + chunk_bytes = '128 MiB' + # chunk_bytes = '64 MiB' + + with dask.config.set({"array.chunk-size": chunk_bytes}) : + # with dask.config.set(**{'array.slicing.split_large_chunks': True}) : + f_obj = open_file_or_filename(self.filename) + self.nc = xr.open_dataset(f_obj, + decode_cf=True, + mask_and_scale=False, + chunks={"xc": "auto", "yc": "auto"}) + + self.metadata = {} + + self.metadata["start_time"] = self.get_start_time() + self.metadata["end_time"] = self.get_end_time() + + print("Reading: {}".format(filename), " start: {} ".format(self.start_time), \ + " end: {}".format(self.end_time)) + self._cache = {} + + self.sensor_name() + self.res() + self.longitudeReelle = self.satlon() + self.longitudedeprojection = self.projlon() + + self.zone = self.nc.attrs["Area_of_acquisition"] + # globe, europe. + + # Reading the needed attributes. + self.initialisation_dataset() + # __init__() + + + def _chunk_bytes_for_resolution(self) -> int: + """Get a best-guess optimal chunk size for resolution-based chunking. + + First a chunk size is chosen for the provided Dask setting `array.chunk-size` + and then aligned with a hardcoded on-disk chunk size of 226. This is then + adjusted to match the current resolution. + + This should result in 500 meter data having 4 times as many pixels per + dask array chunk (2 in each dimension) as 1km data and 8 times as many + as 2km data. As data is combined or upsampled geographically the arrays + should not need to be rechunked. Care is taken to make sure that array + chunks are aligned with on-disk file chunks at all resolutions, but at + the cost of flexibility due to a hardcoded on-disk chunk size of 226 + elements per dimension. + + """ + num_high_res_elems_per_dim = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats + # assume on-disk chunk size of 226 + # this is true for all CSPP Geo GRB output (226 for all sectors) and full disk from other sources + # 250 has been seen for AWS/CLASS CONUS, Mesoscale 1, and Mesoscale 2 files + # we align this with 4 on-disk chunks at 500m, so it will be 2 on-disk chunks for 1km, and 1 for 2km + high_res_elems_disk_aligned = round(max(num_high_res_elems_per_dim / (4 * 226), 1)) * (4 * 226) + low_res_factor = int(self.filetype_info.get("resolution", 2000) // 500) + res_elems_per_dim = int(high_res_elems_disk_aligned / low_res_factor) + + print("num_high_res_elems_per_dim = ", num_high_res_elems_per_dim, \ + " res_elems_per_dim = ", res_elems_per_dim) + # num_high_res_elems_per_dim = 5792.618751480198 res_elems_per_dim = 1356 + return (res_elems_per_dim ** 2) * 2 # 16-bit integers on disk + + + def sensor_name(self) : + """Get the sensor name.""" + # the sensor and platform names are stored together, eg: MSG1/SEVIRI + variable = self.nc["satellite"] + self.plateforme = variable.attrs["id"] + + if "goes" in self.plateforme : + self.sensor = "abi" + elif "msg" in self.plateforme : + self.sensor = "seviri" + elif "MSG" in self.plateforme : + self.sensor = "seviri" + elif "mtg" in self.plateforme : + self.sensor = "fci" + elif "hima" in self.plateforme : + self.sensor = "ahi" + elif "MTG" in self.plateforme : + self.sensor = "fci" + else : + raise NameError("Unsupported satellite platform : " + self.plateforme) + + # Icare and météo france use non-standard platform names. Change is needed for pyspectral. + # pyspectral/rsr_seviri_Meteosat-10.h5 and not rsr_seviri_msg3_seviri.h5 in the call + # Calculator(metadata["platform_name"], metadata["sensor"], metadata["name"]). + if self.plateforme == "msg1" or self.plateforme == "msg01" or self.plateforme == "MSG1" : + self.plateforme = "Meteosat-08" + elif self.plateforme == "msg2" or self.plateforme == "msg02" or self.plateforme == "MSG2" : + self.plateforme = "Meteosat-09" + elif self.plateforme == "msg3" or self.plateforme == "msg03" or self.plateforme == "MSG3" : + self.plateforme = "Meteosat-10" + elif self.plateforme == "msg4" or self.plateforme == "msg04" or self.plateforme == "MSG4" : + self.plateforme = "Meteosat-11" + elif self.plateforme == "mtgi1" or self.plateforme == "mtg1" or self.plateforme == "MTG01" : + self.plateforme = "Meteosat-12" + elif self.plateforme == "goes16" : + self.plateforme = "GOES-16" + elif self.plateforme == "goes17" : + self.plateforme = "GOES-17" + elif self.plateforme == "goes18" : + self.plateforme = "GOES-18" + elif self.plateforme == "goes19" : + self.plateforme = "GOES-19" + elif self.plateforme == "hima08" : + self.plateforme = "Himawari-8" + elif self.plateforme == "hima09" : + self.plateforme = "Himawari-9" + # sensor_name() + + + + def satlon(self): + """Get the satellite longitude.""" + variable = self.nc["satellite"] + longitudeReelle = variable.attrs["lon"] + return longitudeReelle + + + def projlon(self): + """Get the projection longitude.""" + variable = self.nc["geos"] + longitudedeprojection = variable.attrs["longitude_of_projection_origin"] + return longitudedeprojection + + + @property + def projection(self): + """Get the projection.""" + return "geos" + + + def res(self) : + """Get the resolution.""" + # The resolution can be read in the attribute geotransform of the follonwing variables : + # GeosCoordinateSystem500m, GeosCoordinateSystem_h, + # GeosCoordinateSystem1km, GeosCoordinateSystem2km, + # GeosCoordinateSystem. + # cfac, lfac, coff, loff can be read in the variables ImageNavigationxxx. + trouve = False + try : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem500m"] + variableX = self.nc["X500m"] + variableY = self.nc["Y500m"] + chaineNavigation = "ImageNavigation500m" + trouve = True + except KeyError : + None + if not trouve : + try : + # Hrv from msg. + variable = self.nc["GeosCoordinateSystem_h"] + variableX = self.nc["X_h"] + variableY = self.nc["Y_h"] + chaineNavigation = "ImageNavigation_h" + trouve = True + except KeyError : + None + if not trouve : + try : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem1km"] + variableX = self.nc["X1km"] + variableY = self.nc["Y1km"] + chaineNavigation = "ImageNavigation1km" + trouve = True + except KeyError : + None + if not trouve : + try : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem2km"] + variableX = self.nc["X2km"] + variableY = self.nc["Y2km"] + chaineNavigation = "ImageNavigation2km" + trouve = True + except KeyError : + None + + if not trouve : + try : + # Msg in 3 kms. + variable = self.nc["GeosCoordinateSystem"] + variableX = self.nc["X"] + variableY = self.nc["Y"] + chaineNavigation = "ImageNavigation" + trouve = True + except KeyError : + None + + if not trouve : + print("Variables GeosCoordinateSystemXX not founded.") + exit(1) + + geotransform = variable.attrs["GeoTransform"] + + # print("geotransform = ", geotransform) + # geotransform = -5570254, 3000.40604, 0, 5570254, 0, -3000.40604 + morceaux = geotransform.split(", ") + self.resolution = float(morceaux[1]) # 3000.40604 + # print("resolution = ", self.resolution) + + self.X = variableX[:] + self.nbpix = self.X.shape[0] + self.Y = variableY[:] + self.nblig = self.Y.shape[0] + + variable = self.nc[chaineNavigation] + self.cfac = float(variable.attrs["CFAC"]) + self.lfac = float(variable.attrs["LFAC"]) + self.coff = float(variable.attrs["COFF"]) + self.loff = float(variable.attrs["LOFF"]) + # res() + + + def get_end_time(self) : + """Get the end time. Global attribute of the netcdf.""" + attr = self.nc.attrs["time_coverage_end"] + # YYYY-MM-DDTHH:MM:SSZNNN + # In some versions milliseconds are present, sometimes not. + # For the goesr : 2024-06-28T10:00:21.1Z + longueurchaine = len(attr) + if longueurchaine == 22 : + # Goesr. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") + elif longueurchaine == 20 : + # Mtg. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") + else : + # Msg, hima. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ%f") + return stacq + # get_end_time() + + def get_start_time(self) : + """Get the start time. Global attribute of the netcdf.""" + attr = self.nc.attrs["time_coverage_start"] + + # YYYY-MM-DDTHH:MM:SSZ + # In some versions milliseconds are present, sometimes not. + longueurchaine = len(attr) + if longueurchaine == 22 : + # Goesr. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") + elif longueurchaine == 20 : + # Mtg. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") + else : + # Msg, hima. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ%f") + return stacq + # get_start_time() + + @property + def start_time(self) : + return(self.metadata["start_time"]) + @property + def end_time(self) : + return(self.metadata["end_time"]) + + + @property + def alt(self): + """Get the altitude.""" + variable = self.nc["satellite"] + altitude = variable.attrs["dst"] + # 36000000 meters. + altitude += 6378169. # equatorial radius of the earth. + return altitude + + + # ----------------------------------------------------- # + # ----------------------------------------------------- # + def preparer_metadata(self, variable) : + """Get the metadata.""" + mda = {} + + attributs = variable.attrs + for name in attributs : + mda.update({name: attributs.get(name)}) + + mda.update({ + "start_time": self.start_time, + "end_time": self.end_time, + "platform_name": self.plateforme, + "sensor": self.sensor, + "zone": self.zone, + "projection_altitude": self.alt, + "cfac": self.cfac, + "lfac": self.lfac, + "coff": self.coff, + "loff": self.loff, + "resolution": self.resolution, + "satellite_actual_longitude": self.longitudeReelle, + "projection_longitude": self.longitudedeprojection, + "projection_type": self.projection + }) + + # Placer ici les paramètres orbitaux, une seule fois ? + mda.update(self.orbital_param()) + + return mda + # preparer_metadata(). + + + def _get_dsname(self, ds_id) : + """Return the correct dataset name based on requested band.""" + # ds_id = DataID(name='vis_08', wavelength=WavelengthRange(min=0.815, central=0.865, max=0.915, unit='µm'), + # resolution=2000, calibration=, modifiers=()) + + # For mtg : + if ds_id["name"] == 'vis_04' : # Name in satpy. + ds_get_name = "VIS004" # Name in icare/meteofrance netcdf. + elif ds_id["name"] == 'vis_05' : + ds_get_name = "VIS005" + elif ds_id["name"] == 'vis_06' : + ds_get_name = "VIS006" + elif ds_id["name"] == 'vis_08' : + ds_get_name = "VIS008" + elif ds_id["name"] == 'vis_09' : + ds_get_name = "VIS009" + elif ds_id["name"] == 'nir_13' : + ds_get_name = "IR_013" + elif ds_id["name"] == 'nir_16' : + ds_get_name = "IR_016" + elif ds_id["name"] == 'nir_22' : + ds_get_name = "IR_022" + elif ds_id["name"] == 'ir_38' : + ds_get_name = "IR_038" + elif ds_id["name"] == 'wv_63' : + ds_get_name = "WV_063" + elif ds_id["name"] == 'wv_73' : + ds_get_name = "WV_073" + elif ds_id["name"] == 'ir_87' : + ds_get_name = "IR_087" + elif ds_id["name"] == 'ir_97' : + ds_get_name = "IR_097" + elif ds_id["name"] == 'ir_105' : + ds_get_name = "IR_105" + elif ds_id["name"] == 'ir_123' : + ds_get_name = "IR_123" + elif ds_id["name"] == 'ir_133' : + ds_get_name = "IR_133" + + # For msg, the satpy and icare channel names are the same. + elif ds_id["name"] == 'IR_039' or ds_id["name"] == 'IR_016' or ds_id["name"] == 'VIS008' \ + or ds_id["name"] == 'IR_087' or ds_id["name"] == 'IR_097' or ds_id["name"] == 'IR_108' \ + or ds_id["name"] == 'IR_120' or ds_id["name"] == 'IR_134' or ds_id["name"] == 'VIS006' \ + or ds_id["name"] == 'WV_062' or ds_id["name"] == 'WV_073' or ds_id["name"] == 'HRV' : + ds_get_name = ds_id["name"] + + # For the goesr : + elif ds_id["name"] == 'C01' : + ds_get_name = "VIS_004" + elif ds_id["name"] == 'C02' : + ds_get_name = "VIS_006" + elif ds_id["name"] == 'C03' : + ds_get_name = "VIS_008" + elif ds_id["name"] == 'C04' : + ds_get_name = "VIS_014" + elif ds_id["name"] == 'C05' : + ds_get_name = "VIS_016" + elif ds_id["name"] == 'C06' : + ds_get_name = "VIS_022" + elif ds_id["name"] == 'C07' : + ds_get_name = "IR_039" + elif ds_id["name"] == 'C08' : + ds_get_name = "IR_062" + elif ds_id["name"] == 'C09' : + ds_get_name = "IR_069" + elif ds_id["name"] == 'C10' : + ds_get_name = "IR_073" + elif ds_id["name"] == 'C11' : + ds_get_name = "IR_085" + elif ds_id["name"] == 'C12' : + ds_get_name = "IR_096" + elif ds_id["name"] == 'C13' : + ds_get_name = "IR_103" + elif ds_id["name"] == 'C14' : + ds_get_name = "IR_114" + elif ds_id["name"] == 'C15' : + ds_get_name = "IR_123" + elif ds_id["name"] == 'C16' : + ds_get_name = "IR_133" + + # For himawari. + elif ds_id["name"] == 'B01' : # Name in satpy. + ds_get_name = "VIS004" # Name in icare/meteofrance netcdf. + elif ds_id["name"] == 'B02' : + ds_get_name = "VIS005" + elif ds_id["name"] == 'B03' : + ds_get_name = "VIS006" + elif ds_id["name"] == 'B04' : + ds_get_name = "VIS008" + elif ds_id["name"] == 'B05' : + ds_get_name = "IR_016" + elif ds_id["name"] == 'B06' : + ds_get_name = "IR_022" + elif ds_id["name"] == 'B07' : + ds_get_name = "IR_038" + elif ds_id["name"] == 'B08' : + ds_get_name = "WV_062" + elif ds_id["name"] == 'B09' : + ds_get_name = "WV_069" + elif ds_id["name"] == 'B10' : + ds_get_name = "WV_073" + elif ds_id["name"] == 'B11' : + ds_get_name = "IR_085" + elif ds_id["name"] == 'B12' : + ds_get_name = "IR_096" + elif ds_id["name"] == 'B13' : + ds_get_name = "IR_104" + elif ds_id["name"] == 'B14' : + ds_get_name = "IR_112" + elif ds_id["name"] == 'B15' : + ds_get_name = "IR_123" + elif ds_id["name"] == 'B16' : + ds_get_name = "IR_132" + + else : + print("Soft not adaptated for this channel : ds_id = ", ds_id["name"]) + exit(1) + + return ds_get_name + # _get_dsname() + + + def initialisation_dataset(self) : + listeToutesVariables = {"VIS004", "VIS_004", "VIS005", "VIS_005", "VIS006", "VIS_006", "VIS006", \ + "VIS_008", "VIS008", "VIS_009", "VIS009", "IR_013", "IR_016", "IR_022", "IR_038", "IR_039", \ + "IR_062", "WV_062", "WV_063", "IR_069", "WV_069", "IR_073", "WV_073", "IR_085", "IR_087", \ + "IR_096", "IR_097", "IR_103", "IR_104", "IR_105", "IR_108", "IR_112", "IR_120", "IR_123", \ + "IR_132", "IR_133", "IR_134", "HRV"} + + self.mda = {} + self.scale_factor = {} + self.offset = {} + self.alpha = {} + self.beta = {} + self.bandfactor = {} + self.nuc = {} + self.variableComptes = {} + + # Loop over the all possible ds_get_name of every satellites : + for ds_get_name in listeToutesVariables : + if ds_get_name in self.nc : + + variable = self.nc[ds_get_name] + attributs = variable.attrs + + self.scale_factor[ds_get_name] = attributs["scale_factor"] + self.offset[ds_get_name] = attributs["add_offset"] + + if "nuc" in attributs : + # Brightness temperature. + self.alpha[ds_get_name] = attributs["alpha"] + self.beta[ds_get_name] = attributs["beta"] + + nomRetourComptes = "Temp_to_Native_count_" + ds_get_name + + elif "bandfactor" in attributs : + # Albedo. + self.bandfactor[ds_get_name] = attributs["bandfactor"] + + nomRetourComptes = "Albedo_to_Native_count_" + ds_get_name + + else : + print("Nuc or bandfactor not founded int the attributs of ", ds_get_name) + exit(1) + + self.variableComptes[ds_get_name] = self.nc[nomRetourComptes] + # (65536). Correspondence from 0 to 65535 towards the original spatial agency counts. + + self.mda[ds_get_name] = self.preparer_metadata(variable) + # initialisation_dataset() + + + def get_dataset(self, ds_id, ds_info) : + """Get the dataset. + ds_id["calibration"] = key["calibration"] + = ["brightness_temperature", "reflectance", "radiance", "counts"] + """ + ds_get_name = self._get_dsname(ds_id) # "IR_096" + + variable = self.nc[ds_get_name] + + scale_factor = self.scale_factor[ds_get_name] + offset = self.offset[ds_get_name] + + # print(ds_get_name, ds_id["calibration"], " min = ", np.min(variable[:].values), + # " max = ", np.max(variable[:].values), "dims = ", variable.dims) + # WV_062 calibration.brightness_temperature, from -9000 to 4000 + # VIS006 calibration.reflectance, from 0 to 10000 + + mda = self.mda[ds_get_name] + mda.update(ds_info) + + calibration = ds_id["calibration"] + if calibration == "counts" : + # Come back to the original counts of the hrit... + + # Variable is between -9000 to 4000 (temperature) + # or between 0 to 10000 (albedo). + + variable += 32768 # 0 to 65535 + + if offset == 0. : + # Albedo. + nom = "Albedo_to_Native_count_" + ds_get_name + else : + nom = "Temp_to_Native_count_" + ds_get_name + + variableComptes = self.variableComptes[ds_get_name] + # (65536). Correspondence from 0 to 65535 towards the original spatial agency counts. + + # Come back to the original counts of the hrit... + sortie = variablesComptes[AlbedoBT] + + elif calibration == "radiance" : + # Come back to the radiance. + + if offset == 0. : + # Visible channel. + bandfactor = self.bandfactor[ds_get_name] + + # Variable is an albedo from 0 to 10000. + variable = variable * scale_factor/ 100. * bandfactor + # => variable is a reflectance between 0 and 1. + # radiance in mWm-2sr-1(cm-1)-1 + else : + # Brightness temperature. + nuc = self.nuc[ds_get_name] + alpha = self.alpha[ds_get_name] + beta = self.beta[ds_get_name] + + variable = variable * scale_factor + offset + # variable is in Kelvin. + variable = variable * alpha + beta + + c1 = 1.1910427e-5 # Planck + c2 = 1.4387752 # Planck + resul1 = c1 * np.power(nuc, 3.) + resul2 = c2 * nuc + val2 = np.exp(resul2 / variable) - 1. + variable = resul1 / val2 + # Radiance in mWm-2sr-1(cm-1)-1 + + elif calibration == "brightness_temperature" : + if offset != 273.15 : + print(ds_get_name, ". offset = ", offset) + print("Soft not intended for a reflectance with a wave length more than 3 microns.") + exit(1) + + variable = variable * scale_factor + offset + # variable is in Kelvin. + + elif calibration == "reflectance" : + if offset != 0. : + print(ds_get_name, ". offset = ", offset) + message = "Soft not intended for a brightness temperature " + message += "with a wave length less than 3 microns." + print(message) + exit(1) + + variable = variable * scale_factor + # variable is an albedo between 0 and 100. + + else : + print("Calibration mode not expected : ", calibration) + exit(1) + + variable = variable.rename({variable.dims[0] : "y", variable.dims[1] : "x"}) + variable.attrs.update(mda) + return variable + # get_dataset() + + + def orbital_param(self) : + orb_param_dict = { + "orbital_parameters": { + "satellite_actual_longitude": self.longitudeReelle, + "satellite_actual_latitude": 0., + "satellite_actual_altitude": self.alt, + "satellite_nominal_longitude": self.longitudedeprojection, + "satellite_nominal_latitude": 0, + "satellite_nominal_altitude": self.alt, + "projection_longitude": self.longitudedeprojection, + "projection_latitude": 0., + "projection_altitude": self.alt, + }} + + return orb_param_dict + + def get_area_def(self, ds_id) : + """Get the area def.""" + + pdict = {} + pdict["cfac"] = np.int32(self.cfac) + pdict["lfac"] = np.int32(self.lfac) + pdict["coff"] = np.float32(self.coff) + pdict["loff"] = np.float32(self.loff) + + pdict["a"] = 6378169 + pdict["b"] = 6356583.8 + pdict["h"] = self.alt - pdict["a"] # 36000000 mètres. + pdict["ssp_lon"] = self.longitudedeprojection + pdict["ncols"] = self.nblig + pdict["nlines"] = self.nbpix + pdict["sweep"] = "y" + + # Force scandir to SEVIRI default, not known from file + pdict["scandir"] = "S2N" + pdict["a_name"] = "geosmsg" + + if self.sensor == "seviri" : + # msg. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geosmsg" + if self.nbpix == 3712 : + pdict["a_desc"] = "MSG/SEVIRI low resolution channel area" + pdict["p_id"] = "msg_lowres" + elif self.nbpix == 11136 : + pdict["a_desc"] = "MSG/SEVIRI HRV channel area" + pdict["p_id"] = "msg_hires" + else : + print("ERROR : not expected resolution for msg : ", self.nbpix) + exit(1) + + elif self.sensor == "fci" : + # mtg. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geosmtg" + if self.nbpix == 5568 : + pdict["a_desc"] = "MTG 2km channel area" + pdict["p_id"] = "mtg_lowres" + elif self.nbpix == 11136 : + pdict["a_desc"] = "MTG 1km channel area" + pdict["p_id"] = "mtg_midres" + elif self.nbpix == 22272 : + pdict["a_desc"] = "MTG 500m channel area" + pdict["p_id"] = "mtg_hires" + else : + print("ERROR : not expected resolution for mtg : ", self.nbpix) + exit(1) + + elif self.sensor == "ahi" : + # Himawari. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geoshima" + if self.nbpix == 5500 : + pdict["a_desc"] = "HIMA 2km channel area" + pdict["p_id"] = "hima_lowres" + elif self.nbpix == 11000 : + pdict["a_desc"] = "HIMA 1km channel area" + pdict["p_id"] = "hima_midres" + elif self.nbpix == 22000 : + pdict["a_desc"] = "HIMA 500m channel area" + pdict["p_id"] = "hima_hires" + else : + print("ERROR : not expected resolution for hima : ", self.nbpix) + exit(1) + + elif self.sensor == "abi" : + # Goesr. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geosgoesr" + pdict["sweep"] = "x" + if self.nbpix == 5424 : + pdict["a_desc"] = "GOESR 2km channel area" + pdict["p_id"] = "goesr_lowres" + elif self.nbpix == 10848 : + pdict["a_desc"] = "GOESR 1km channel area" + pdict["p_id"] = "goesr_midres" + elif self.nbpix == 21696 : + pdict["a_desc"] = "GOESR 500m channel area" + pdict["p_id"] = "goesr_hires" + else : + print("ERROR : not expected resolution for goesr : ", self.nbpix) + exit(1) + else : + print("ERROR : " + self.sensor + " not expected.") + exit(1) + + aex = get_area_extent(pdict) + area = get_area_definition(pdict, aex) + + # print("area = ", area) + return area + # get_area_def() + + + From 2dbc7bea0338e6435919dbd94086721b80a5a2dd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 29 Sep 2024 03:35:39 +0000 Subject: [PATCH 2/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- satpy/etc/readers/goesr_netcdficare.yaml | 3 -- satpy/etc/readers/hima_netcdficare.yaml | 1 - satpy/etc/readers/msg_netcdficare.yaml | 14 ++++---- satpy/etc/readers/mtg_netcdficare.yaml | 2 -- satpy/readers/geos_netcdficare.py | 42 +++++++++++------------- 5 files changed, 25 insertions(+), 37 deletions(-) diff --git a/satpy/etc/readers/goesr_netcdficare.yaml b/satpy/etc/readers/goesr_netcdficare.yaml index 7118e0dd07..4b36914b67 100644 --- a/satpy/etc/readers/goesr_netcdficare.yaml +++ b/satpy/etc/readers/goesr_netcdficare.yaml @@ -304,6 +304,3 @@ datasets: standard_name: counts units: "1" file_type: GOES2km - - - diff --git a/satpy/etc/readers/hima_netcdficare.yaml b/satpy/etc/readers/hima_netcdficare.yaml index cdd3bcfd96..bd3bd70f92 100644 --- a/satpy/etc/readers/hima_netcdficare.yaml +++ b/satpy/etc/readers/hima_netcdficare.yaml @@ -258,4 +258,3 @@ datasets: standard_name: counts units: 1 file_type: [HIMA2km] - diff --git a/satpy/etc/readers/msg_netcdficare.yaml b/satpy/etc/readers/msg_netcdficare.yaml index 01fa4245cb..8885625321 100644 --- a/satpy/etc/readers/msg_netcdficare.yaml +++ b/satpy/etc/readers/msg_netcdficare.yaml @@ -18,10 +18,10 @@ file_types: MSG3km : file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE - file_patterns: ['Mmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + file_patterns: ['Mmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', 'Imultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', 'Mrsmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', - 'Mmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Mmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', 'Imultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', 'Mrsmultic3kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc'] @@ -29,11 +29,11 @@ file_types: MSG1km : file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE - file_patterns: [ 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', - 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + file_patterns: [ 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', - 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', - 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc'] # platform_shortname:5s : msg01, ..., msg04 @@ -170,5 +170,3 @@ datasets: standard_name: brightness_temperature units: "K" file_type: MSG3km - - diff --git a/satpy/etc/readers/mtg_netcdficare.yaml b/satpy/etc/readers/mtg_netcdficare.yaml index 2562feafd7..9b8db1478c 100644 --- a/satpy/etc/readers/mtg_netcdficare.yaml +++ b/satpy/etc/readers/mtg_netcdficare.yaml @@ -331,5 +331,3 @@ datasets: brightness_temperature: standard_name: toa_brightness_temperature units: "K" - - diff --git a/satpy/readers/geos_netcdficare.py b/satpy/readers/geos_netcdficare.py index 31bae413b9..018c84d777 100644 --- a/satpy/readers/geos_netcdficare.py +++ b/satpy/readers/geos_netcdficare.py @@ -100,22 +100,21 @@ import datetime as dt +# from satpy.utils import get_dask_chunk_size_in_bytes +import dask import numpy as np - -from satpy.readers._geos_area import get_area_definition, get_area_extent - +import xarray as xr from netCDF4 import Dataset from xarray import DataArray -import xarray as xr - -from satpy.readers.file_handlers import BaseFileHandler from satpy._compat import cached_property -# from satpy.utils import get_dask_chunk_size_in_bytes -import dask -from satpy.readers import open_file_or_filename +from satpy.readers._geos_area import get_area_definition, get_area_extent +from satpy.readers.file_handlers import BaseFileHandler import math +from satpy.readers import open_file_or_filename + + class NETCDF_ICARE(BaseFileHandler) : # Cf readers/file_handlers.py. @@ -156,7 +155,7 @@ def __init__(self, filename, filename_info, filetype_info) : self.zone = self.nc.attrs["Area_of_acquisition"] # globe, europe. - + # Reading the needed attributes. self.initialisation_dataset() # __init__() @@ -213,7 +212,7 @@ def sensor_name(self) : self.sensor = "fci" else : raise NameError("Unsupported satellite platform : " + self.plateforme) - + # Icare and météo france use non-standard platform names. Change is needed for pyspectral. # pyspectral/rsr_seviri_Meteosat-10.h5 and not rsr_seviri_msg3_seviri.h5 in the call # Calculator(metadata["platform_name"], metadata["sensor"], metadata["name"]). @@ -266,7 +265,7 @@ def projection(self): def res(self) : """Get the resolution.""" # The resolution can be read in the attribute geotransform of the follonwing variables : - # GeosCoordinateSystem500m, GeosCoordinateSystem_h, + # GeosCoordinateSystem500m, GeosCoordinateSystem_h, # GeosCoordinateSystem1km, GeosCoordinateSystem2km, # GeosCoordinateSystem. # cfac, lfac, coff, loff can be read in the variables ImageNavigationxxx. @@ -441,7 +440,7 @@ def _get_dsname(self, ds_id) : """Return the correct dataset name based on requested band.""" # ds_id = DataID(name='vis_08', wavelength=WavelengthRange(min=0.815, central=0.865, max=0.915, unit='µm'), # resolution=2000, calibration=, modifiers=()) - + # For mtg : if ds_id["name"] == 'vis_04' : # Name in satpy. ds_get_name = "VIS004" # Name in icare/meteofrance netcdf. @@ -611,7 +610,7 @@ def initialisation_dataset(self) : def get_dataset(self, ds_id, ds_info) : """Get the dataset. - ds_id["calibration"] = key["calibration"] + ds_id["calibration"] = key["calibration"] = ["brightness_temperature", "reflectance", "radiance", "counts"] """ ds_get_name = self._get_dsname(ds_id) # "IR_096" @@ -621,7 +620,7 @@ def get_dataset(self, ds_id, ds_info) : scale_factor = self.scale_factor[ds_get_name] offset = self.offset[ds_get_name] - # print(ds_get_name, ds_id["calibration"], " min = ", np.min(variable[:].values), + # print(ds_get_name, ds_id["calibration"], " min = ", np.min(variable[:].values), # " max = ", np.max(variable[:].values), "dims = ", variable.dims) # WV_062 calibration.brightness_temperature, from -9000 to 4000 # VIS006 calibration.reflectance, from 0 to 10000 @@ -633,7 +632,7 @@ def get_dataset(self, ds_id, ds_info) : if calibration == "counts" : # Come back to the original counts of the hrit... - # Variable is between -9000 to 4000 (temperature) + # Variable is between -9000 to 4000 (temperature) # or between 0 to 10000 (albedo). variable += 32768 # 0 to 65535 @@ -653,13 +652,13 @@ def get_dataset(self, ds_id, ds_info) : elif calibration == "radiance" : # Come back to the radiance. - if offset == 0. : + if offset == 0. : # Visible channel. bandfactor = self.bandfactor[ds_get_name] # Variable is an albedo from 0 to 10000. variable = variable * scale_factor/ 100. * bandfactor - # => variable is a reflectance between 0 and 1. + # => variable is a reflectance between 0 and 1. # radiance in mWm-2sr-1(cm-1)-1 else : # Brightness temperature. @@ -745,7 +744,7 @@ def get_area_def(self, ds_id) : # Force scandir to SEVIRI default, not known from file pdict["scandir"] = "S2N" pdict["a_name"] = "geosmsg" - + if self.sensor == "seviri" : # msg. pdict["scandir"] = "N2S" @@ -776,7 +775,7 @@ def get_area_def(self, ds_id) : else : print("ERROR : not expected resolution for mtg : ", self.nbpix) exit(1) - + elif self.sensor == "ahi" : # Himawari. pdict["scandir"] = "N2S" @@ -821,6 +820,3 @@ def get_area_def(self, ds_id) : # print("area = ", area) return area # get_area_def() - - - From 93987374864b41753af2c0ae887500ad86a202b3 Mon Sep 17 00:00:00 2001 From: lperier Date: Mon, 14 Oct 2024 16:52:43 +0200 Subject: [PATCH 3/9] Add files via upload Pytest program to test the reader geos_netcdficare.py. --- test_geos_netcdficare.py | 308 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 308 insertions(+) create mode 100644 test_geos_netcdficare.py diff --git a/test_geos_netcdficare.py b/test_geos_netcdficare.py new file mode 100644 index 0000000000..cf1f5401ad --- /dev/null +++ b/test_geos_netcdficare.py @@ -0,0 +1,308 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2019 Satpy developers +# +# This file is part of satpy. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with +# satpy. If not, see . + +"""Tests for the Icare MeteoFrance netcdfs reader, + satpy/readers/geos_netcdficare.py. + +SYNTAXE : + pytest + + Before, copy the satpy/tests/reader_tests/test_geos_netcdfcare.py file + into the pytest directory pytest_projec. + +EXECUTION TIME : + 20 seconds. +DATE OF CREATION : + 2024 11th october. +LAST VERSIONS : + +AUTHOR : + Meteo France. +""" + +import os + +import numpy as np + +from satpy.scene import Scene +from satpy import find_files_and_readers + +from datetime import datetime +from netCDF4 import Dataset + + +class TestGeosNetcdfIcareReader() : + # Test of the geos_netcdficare reader. + # This reader has been build for the Icare Meteo France netcdfs. + + def test_geos_netcdficare(self, tmp_path) : + """ A dummy netcdf is built. + A scene self.scn for the convection product for the same date + is built. We check that the scene parameters are the same + as thoses in the dummy netcdf. + This procedure is called by pytest. + """ + + self.init(tmp_path) + + startTime = self.scn.start_time + startTimeString = startTime.strftime('%Y-%m-%dT%H:%M:%S') + # 2024-06-28T10:00:40 + assert startTimeString == self.expectedStartTime + + endTime = self.scn.end_time + endTimeString = endTime.strftime('%Y-%m-%dT%H:%M:%S') + # 2024-06-28T10:12:41 + assert endTimeString == self.expectedEndTime + + sensor = self.scn.sensor_names + for isensor in sensor : + capteur = isensor + assert capteur == self.expectedSensor + + platform = "erreur" + altitude = -1. + longitude = 999. + + for data_arr in self.values : + # values come from the scene. + # print("data_arr.attrs = ", data_arr.attrs) + if "platform_name" in data_arr.attrs : + platform = data_arr.attrs["platform_name"] + if "orbital_parameters" in data_arr.attrs : + subAttr = data_arr.attrs["orbital_parameters"] + if "satellite_actual_altitude" in subAttr : + altitude = subAttr["satellite_actual_altitude"] + if "satellite_actual_longitude" in data_arr.attrs : + longitude = data_arr.attrs["satellite_actual_longitude"] + + longitude = float(int(longitude * 10.)) / 10. + assert platform == self.expectedPlatform + assert longitude == self.expectedLongitude + assert altitude == self.expectedAltitude + + xr = self.scn.to_xarray_dataset() + # print("xr = ", xr) + matrice = xr["convection"] + nblin = matrice.shape[1] + nbpix = matrice.shape[2] + assert nblin == self.expectedNblin + assert nbpix == self.expectedNbpix + + cfac = xr.attrs["cfac"] + assert cfac == self.expectedCfac + lfac = xr.attrs["lfac"] + assert lfac == self.expectedLfac + coff = xr.attrs["coff"] + assert coff == self.expectedCoff + loff = xr.attrs["loff"] + assert loff == self.expectedLoff + + satpyId = xr.attrs["_satpy_id"] + # DataID(name='convection', resolution=3000.403165817) + # Cf satpy/dataset/dataid.py. + + resolution = satpyId.get("resolution") + resolution = float(int(resolution * 10.)) / 10. + assert resolution == self.expectedResolution + + # A picture of convection composite will be displayed. + self.scn.show("convection") + print("The picture should be pink.") + # test_geos_netcdficare(self, tmp_path) + + def init(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_geos_netcdficare(). + """ + self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" + self.filepath = tmp_path + + self.buildNetcdf(self.netcdfName) + + # We will check that the parameters written in the dummy netcdf can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:12:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "Meteosat-10" + self.expectedSensor = "seviri" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 1.3642337E7 + self.expectedLfac = 1.3642337E7 + self.expectedCoff = 1857.0 + self.expectedLoff = 1857.0 + self.expectedResolution = 3000.4 + self.expectedNblin = 3712 + self.expectedNbpix = 3712 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for in the filepath directory. + yaml_file = 'msg_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + print("Found myfiles = ", myfiles) + # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + + print(self.scn.available_dataset_names()) + # ['IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', + # 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073'] + + print(self.scn.available_composite_names()) + """ Static compositor {'_satpy_id': DataID(name='_night_background'), + 'standard_name': 'night_background', + 'prerequisites': [], + 'optional_prerequisites': []} + Static compositor {'_satpy_id': DataID(name='_night_background_hires'), + 'standard_name': 'night_background', 'prerequisites': [], + 'optional_prerequisites': []} + ['airmass', 'ash', 'cloud_phase_distinction', 'cloud_phase_distinction_raw', + 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection'... + """ + + self.scn.load(['convection']) + self.values = self.scn.values() + # init() + + def buildNetcdf(self, ncName) : + """ + ncName : tmp_path / Mmultic3kmNC4_msg03_202406281000.nc + A dummy icare Meteo France netcdf is built here. + Called by init(). + """ + if os.path.exists(ncName) : + os.remove(ncName) + ncfileOut = Dataset( + ncName, mode="w", clobber=True, + diskless=False, persist=False, format='NETCDF4') + + ncfileOut.createDimension(u'ny', 3712) + ncfileOut.createDimension(u'nx', 3712) + ncfileOut.createDimension(u'numerical_count', 65536) + ncfileOut.setncattr("time_coverage_start", "2024-06-28T10:00:09Z383") + ncfileOut.setncattr("time_coverage_end", "2024-06-28T10:12:41Z365") + ncfileOut.setncattr("Area_of_acquisition", "globe") + + fill_value = -32768 + var = ncfileOut.createVariable( + "satellite", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + + var.setncattr("id", "msg03") + var.setncattr("dst", 35786691.) + var.setncattr("lon", float(0.1)) + + var = ncfileOut.createVariable( + "geos", "c", zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None) + var.setncattr("longitude_of_projection_origin", 0.) + + var = ncfileOut.createVariable( + "GeosCoordinateSystem", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + var.setncattr( + "GeoTransform", + "-5570254, 3000.40604, 0, 5570254, 0, -3000.40604") + + var = ncfileOut.createVariable( + "ImageNavigation", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + var.setncattr("CFAC", 1.3642337E7) + var.setncattr("LFAC", 1.3642337E7) + var.setncattr("COFF", 1857.0) + var.setncattr("LOFF", 1857.0) + + var = ncfileOut.createVariable( + "X", 'float32', u'nx', zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + x0 = -5570254. + dx = 3000.40604 + var[:] = np.array(([(x0 + dx * i) for i in range(3712)])) + + var = ncfileOut.createVariable( + "Y", 'float32', u'ny', zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + y0 = 5570254. + dy = -3000.40604 + var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) + + for channel in {"VIS006", "VIS008", "IR_016"} : + var = ncfileOut.createVariable( + channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', + least_significant_digit=None, fill_value=fill_value) + var[:] = np.array(([[i * 2 for i in range(3712)] for j in range(3712)])) + # Hundredths of albedo between 0 and 10000. + var.setncattr("scale_factor", 0.01) + var.setncattr("add_offset", 0.) + var.setncattr("bandfactor", 20.76) + var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + + var = ncfileOut.createVariable( + "Albedo_to_Native_count_" + channel, 'short', + 'numerical_count', zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None, + fill_value=-9999) + var[:] = np.array(([-9999 for i in range(65536)])) + # In order to come back to the native datas on 10, 12 or 16 bits. + + for channel in { + "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", + "IR_108", "IR_120", "IR_134"} : + var = ncfileOut.createVariable( + channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', + least_significant_digit=None, fill_value=fill_value) + var[:] = np.array( + ([[-9000 + j * 4 for i in range(3712)] for j in range(3712)])) + # Hundredths of celcius degrees. + var.setncattr("scale_factor", 0.01) + var.setncattr("add_offset", 273.15) + var.setncattr("nuc", 1600.548) + var.setncattr("alpha", 0.9963) + var.setncattr("beta", 2.185) + var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + + var = ncfileOut.createVariable( + "Temp_to_Native_count_" + channel, 'short', + 'numerical_count', zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None, + fill_value=-9999) + var[:] = np.array(([-9999 for i in range(65536)])) + # In order to come back to the native datas on 10, 12 or 16 bits. + ncfileOut.close + # buildNetcdf() + + # class TestGeosNetcdfIcareReader. From 551fa37afc38feea47bb4e707c2fc551578a5ed9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?laurent=20p=C3=A9rier?= Date: Mon, 14 Oct 2024 17:42:41 +0200 Subject: [PATCH 4/9] Answer to codescene request. Pytest program for this reader. --- satpy/readers/geos_netcdficare.py | 881 +++++++++--------- .../reader_tests/test_geos_netcdficare.py | 308 ++++++ 2 files changed, 725 insertions(+), 464 deletions(-) create mode 100644 satpy/tests/reader_tests/test_geos_netcdficare.py diff --git a/satpy/readers/geos_netcdficare.py b/satpy/readers/geos_netcdficare.py index 018c84d777..d4cd1cff12 100644 --- a/satpy/readers/geos_netcdficare.py +++ b/satpy/readers/geos_netcdficare.py @@ -21,8 +21,10 @@ Introduction ------------ -The ``geos_netcdficare`` reader reads some geostationnary netcdf build by Meteo France and stored at Icare. -The brightness tempeture and albedo are calibrated. +The ``geos_netcdficare`` reader reads some geostationnary netcdf build by +Meteo France and stored at Icare. + +The brightness tempera ture and albedo are calibrated. That has been stored by the ICARE Data and Services Center Data can be accessed via: http://www.icare.univ-lille1.fr @@ -30,33 +32,41 @@ This reader concerns the following netcdfs : . msg with a longitude near 0° : -Mmultic3kmNC4_msg03_20231113_111500.nc Mmultic1kmNC4_msg03_20231113_111500.nc +Mmultic3kmNC4_msg03_20231113_111500.nc +Mmultic1kmNC4_msg03_20231113_111500.nc . Msg rapid scan with a longitude near 9.5° : -Mrsmultic3kmNC4_msg03_20231113_111500.nc Mrsmultic1kmNC4_msg03_20231113_111500.nc +Mrsmultic3kmNC4_msg03_20231113_111500.nc +Mrsmultic1kmNC4_msg03_20231113_111500.nc . Msg with a longitude near 42° : -Imultic3kmNC4_msg03_20231113_111500.nc Imultic1kmNC4_msg03_20231113_111500.nc +Imultic3kmNC4_msg03_20231113_111500.nc +Imultic1kmNC4_msg03_20231113_111500.nc . Himawari : -Jmultic2kmNC4_hima09_20231113_111500.nc Jmultic1kmNC4_hima09_20231113_111500.nc Jmultic500mNC4_hima09_20231113_111500.nc +Jmultic2kmNC4_hima09_20231113_111500.nc +Jmultic1kmNC4_hima09_20231113_111500.nc +Jmultic500mNC4_hima09_20231113_111500.nc . Goesr near -137° : -Wmultic2kmNC4_goes16_202406281000.nc. The better resolution are not built at Lannion, only at Tahiti. +Wmultic2kmNC4_goes16_202406281000.nc. +The better resolution are not built at Lannion, only at Tahiti. . Goesr in -75° : -Emultic2kmNC4_goes16_202406281000.nc Emultic1kmNC4_goes16_202406281000.nc Emultic500mNC4_goes16_202406281000.nc +Emultic2kmNC4_goes16_202406281000.nc +Emultic1kmNC4_goes16_202406281000.nc +Emultic500mNC4_goes16_202406281000.nc . Mtg : -Mmultic2km_mtgi1_20240104_090000.nc Mmultic1km_mtgi1_20240104_090000.nc Mmultic500m_mtgi1_20240104_090000.nc +Mmultic2km_mtgi1_20240104_090000.nc +Mmultic1km_mtgi1_20240104_090000.nc +Mmultic500m_mtgi1_20240104_090000.nc Example: -------- Here is an example how to read the data in satpy: -.. code-block:: python - from satpy import Scene import glob @@ -64,55 +74,36 @@ scn = Scene(filenames = filenames, reader = 'hima_netcdficare') scn.load(['true_color']) # scn.load(['VIS006']) - my_area = AreaDefinition('my_area', 'zone', 'my_area', - '+proj=latlong +lon_0=0 +a=6378169 +b=6356583 +h=35785831 +x_0=0 +y_0=0 +pm=0', - 8500, 4000, - [-180., -80., 180., 80], - nprocs=16) + my_area = AreaDefinition( + 'my_area', 'zone', 'my_area', + '+proj=latlong +lon_0=0 +a=6378169 +b=6356583 +h=35785831 +x_0=0 + +y_0=0 +pm=0', + 8500, 4000, + [-180., -80., 180., 80], + nprocs=16) natscn = scn.resample(my_area, resampler='nearest') natscn.save_dataset(composite_name, filename = filename_image_out) -Output: - -.. code-block:: none - - - dask.array - Coordinates: - crs object +proj=geos +a=6378169.0 +b=6356583.8 +lon_0=0.0 +h=35785831.0 +units=m +type=crs - * y (y) float64 5.566e+06 5.563e+06 5.56e+06 ... -5.566e+06 -5.569e+06 - * x (x) float64 -5.566e+06 -5.563e+06 -5.56e+06 ... 5.566e+06 5.569e+06 - Attributes: - start_time: 2004-12-29 12:15:00 - end_time: 2004-12-29 12:27:44 - area: Area ID: geosmsg\nDescription: MSG/SEVIRI low resol... - name: IR_108 - resolution: 3000.403165817 - calibration: brightness_temperature - polarization: None - level: None - modifiers: () - ancillary_variables: [] - - """ import datetime as dt -# from satpy.utils import get_dask_chunk_size_in_bytes -import dask import numpy as np -import xarray as xr -from netCDF4 import Dataset -from xarray import DataArray -from satpy._compat import cached_property from satpy.readers._geos_area import get_area_definition, get_area_extent + +# from netCDF4 import Dataset +# from xarray import DataArray +import xarray as xr + from satpy.readers.file_handlers import BaseFileHandler -import math +# from satpy._compat import cached_property +# from satpy.utils import get_dask_chunk_size_in_bytes +import dask from satpy.readers import open_file_or_filename +# import math class NETCDF_ICARE(BaseFileHandler) : @@ -121,30 +112,29 @@ class NETCDF_ICARE(BaseFileHandler) : def __init__(self, filename, filename_info, filetype_info) : """Init the file handler.""" - super().__init__(filename, filename_info, filetype_info) # cache_var_size=0, cache_handle = False) - # self.ncfileIn = Dataset(filename, mode = "r", clobber=True, diskless=False, persist=False, format='NETCDF4') - - # self.filename = filename - # self.filetype_info = filetype_info + super().__init__(filename, filename_info, filetype_info) # chunk_bytes = self._chunk_bytes_for_resolution() chunk_bytes = '128 MiB' # chunk_bytes = '64 MiB' + # with dask.config.set(**{'array.slicing.split_large_chunks': True}) : with dask.config.set({"array.chunk-size": chunk_bytes}) : - # with dask.config.set(**{'array.slicing.split_large_chunks': True}) : f_obj = open_file_or_filename(self.filename) - self.nc = xr.open_dataset(f_obj, - decode_cf=True, - mask_and_scale=False, - chunks={"xc": "auto", "yc": "auto"}) + self.nc = xr.open_dataset( + f_obj, decode_cf=True, mask_and_scale=False, + chunks={"xc": "auto", "yc": "auto"}) self.metadata = {} - self.metadata["start_time"] = self.get_start_time() - self.metadata["end_time"] = self.get_end_time() + self.metadata["start_time"] = self.get_endOrStartTime( + "time_coverage_start") + self.metadata["end_time"] = self.get_endOrStartTime( + "time_coverage_end") - print("Reading: {}".format(filename), " start: {} ".format(self.start_time), \ + print( + "Reading: {}".format(filename), + " start: {} ".format(self.start_time), " end: {}".format(self.end_time)) self._cache = {} @@ -160,168 +150,124 @@ def __init__(self, filename, filename_info, filetype_info) : self.initialisation_dataset() # __init__() - - def _chunk_bytes_for_resolution(self) -> int: - """Get a best-guess optimal chunk size for resolution-based chunking. - - First a chunk size is chosen for the provided Dask setting `array.chunk-size` - and then aligned with a hardcoded on-disk chunk size of 226. This is then - adjusted to match the current resolution. - - This should result in 500 meter data having 4 times as many pixels per - dask array chunk (2 in each dimension) as 1km data and 8 times as many - as 2km data. As data is combined or upsampled geographically the arrays - should not need to be rechunked. Care is taken to make sure that array - chunks are aligned with on-disk file chunks at all resolutions, but at - the cost of flexibility due to a hardcoded on-disk chunk size of 226 - elements per dimension. - - """ - num_high_res_elems_per_dim = math.sqrt(get_dask_chunk_size_in_bytes() / 4) # 32-bit floats - # assume on-disk chunk size of 226 - # this is true for all CSPP Geo GRB output (226 for all sectors) and full disk from other sources - # 250 has been seen for AWS/CLASS CONUS, Mesoscale 1, and Mesoscale 2 files - # we align this with 4 on-disk chunks at 500m, so it will be 2 on-disk chunks for 1km, and 1 for 2km - high_res_elems_disk_aligned = round(max(num_high_res_elems_per_dim / (4 * 226), 1)) * (4 * 226) - low_res_factor = int(self.filetype_info.get("resolution", 2000) // 500) - res_elems_per_dim = int(high_res_elems_disk_aligned / low_res_factor) - - print("num_high_res_elems_per_dim = ", num_high_res_elems_per_dim, \ - " res_elems_per_dim = ", res_elems_per_dim) - # num_high_res_elems_per_dim = 5792.618751480198 res_elems_per_dim = 1356 - return (res_elems_per_dim ** 2) * 2 # 16-bit integers on disk - - def sensor_name(self) : - """Get the sensor name.""" - # the sensor and platform names are stored together, eg: MSG1/SEVIRI + """Get the sensor name. + The sensor and platform names are stored together, eg: MSG1/SEVIRI + """ variable = self.nc["satellite"] self.plateforme = variable.attrs["id"] if "goes" in self.plateforme : self.sensor = "abi" - elif "msg" in self.plateforme : + elif ("msg" in self.plateforme) or ("MSG" in self.plateforme) : self.sensor = "seviri" - elif "MSG" in self.plateforme : - self.sensor = "seviri" - elif "mtg" in self.plateforme : + elif ("mtg" in self.plateforme) or ("MTG" in self.plateforme) : self.sensor = "fci" elif "hima" in self.plateforme : self.sensor = "ahi" - elif "MTG" in self.plateforme : - self.sensor = "fci" else : raise NameError("Unsupported satellite platform : " + self.plateforme) + + # Icare and météo france use non-standard platform names. + # Change is needed for pyspectral : + # pyspectral/rsr_seviri_Meteosat-10.h5 in the call + # Calculator(platform_name, sensor, name). + pdict = {} + pdict["msg1"] = "Meteosat-08" + pdict["msg01"] = "Meteosat-08" + pdict["MSG1"] = "Meteosat-08" + pdict["msg2"] = "Meteosat-09" + pdict["msg02"] = "Meteosat-09" + pdict["MSG2"] = "Meteosat-09" + pdict["msg3"] = "Meteosat-10" + pdict["msg03"] = "Meteosat-10" + pdict["MSG3"] = "Meteosat-10" + pdict["msg4"] = "Meteosat-11" + pdict["msg04"] = "Meteosat-11" + pdict["MSG4"] = "Meteosat-11" + pdict["mtgi1"] = "Meteosat-12" + pdict["mtg1"] = "Meteosat-12" + pdict["MTG01"] = "Meteosat-12" + pdict["goes16"] = "GOES-16" + pdict["goes17"] = "GOES-17" + pdict["goes18"] = "GOES-18" + pdict["goes19"] = "GOES-19" + pdict["hima08"] = "Himawari-8" + pdict["hima09"] = "Himawari-9" + + if self.plateforme in pdict : + self.plateforme = pdict[self.plateforme] + else : + print( + "Unsupported satellite platform : " + + self.plateforme) + exit(1) - # Icare and météo france use non-standard platform names. Change is needed for pyspectral. - # pyspectral/rsr_seviri_Meteosat-10.h5 and not rsr_seviri_msg3_seviri.h5 in the call - # Calculator(metadata["platform_name"], metadata["sensor"], metadata["name"]). - if self.plateforme == "msg1" or self.plateforme == "msg01" or self.plateforme == "MSG1" : - self.plateforme = "Meteosat-08" - elif self.plateforme == "msg2" or self.plateforme == "msg02" or self.plateforme == "MSG2" : - self.plateforme = "Meteosat-09" - elif self.plateforme == "msg3" or self.plateforme == "msg03" or self.plateforme == "MSG3" : - self.plateforme = "Meteosat-10" - elif self.plateforme == "msg4" or self.plateforme == "msg04" or self.plateforme == "MSG4" : - self.plateforme = "Meteosat-11" - elif self.plateforme == "mtgi1" or self.plateforme == "mtg1" or self.plateforme == "MTG01" : - self.plateforme = "Meteosat-12" - elif self.plateforme == "goes16" : - self.plateforme = "GOES-16" - elif self.plateforme == "goes17" : - self.plateforme = "GOES-17" - elif self.plateforme == "goes18" : - self.plateforme = "GOES-18" - elif self.plateforme == "goes19" : - self.plateforme = "GOES-19" - elif self.plateforme == "hima08" : - self.plateforme = "Himawari-8" - elif self.plateforme == "hima09" : - self.plateforme = "Himawari-9" # sensor_name() - - - def satlon(self): + def satlon(self) : """Get the satellite longitude.""" variable = self.nc["satellite"] longitudeReelle = variable.attrs["lon"] return longitudeReelle - def projlon(self): """Get the projection longitude.""" variable = self.nc["geos"] longitudedeprojection = variable.attrs["longitude_of_projection_origin"] return longitudedeprojection - @property def projection(self): """Get the projection.""" return "geos" - def res(self) : - """Get the resolution.""" - # The resolution can be read in the attribute geotransform of the follonwing variables : - # GeosCoordinateSystem500m, GeosCoordinateSystem_h, - # GeosCoordinateSystem1km, GeosCoordinateSystem2km, - # GeosCoordinateSystem. - # cfac, lfac, coff, loff can be read in the variables ImageNavigationxxx. - trouve = False - try : + """Get the resolution. + The resolution can be read in the attribute geotransform + of the following variables : + GeosCoordinateSystem500m, GeosCoordinateSystem_h, + GeosCoordinateSystem1km, GeosCoordinateSystem2km, + GeosCoordinateSystem. + cfac, lfac, coff, loff can be read in the variables ImageNavigationxxx. + """ + + if "GeosCoordinateSystem500m" in self.nc : # Mtg, himawari, goesr. variable = self.nc["GeosCoordinateSystem500m"] variableX = self.nc["X500m"] variableY = self.nc["Y500m"] chaineNavigation = "ImageNavigation500m" - trouve = True - except KeyError : - None - if not trouve : - try : - # Hrv from msg. - variable = self.nc["GeosCoordinateSystem_h"] - variableX = self.nc["X_h"] - variableY = self.nc["Y_h"] - chaineNavigation = "ImageNavigation_h" - trouve = True - except KeyError : - None - if not trouve : - try : - # Mtg, himawari, goesr. - variable = self.nc["GeosCoordinateSystem1km"] - variableX = self.nc["X1km"] - variableY = self.nc["Y1km"] - chaineNavigation = "ImageNavigation1km" - trouve = True - except KeyError : - None - if not trouve : - try : - # Mtg, himawari, goesr. - variable = self.nc["GeosCoordinateSystem2km"] - variableX = self.nc["X2km"] - variableY = self.nc["Y2km"] - chaineNavigation = "ImageNavigation2km" - trouve = True - except KeyError : - None - - if not trouve : - try : - # Msg in 3 kms. - variable = self.nc["GeosCoordinateSystem"] - variableX = self.nc["X"] - variableY = self.nc["Y"] - chaineNavigation = "ImageNavigation" - trouve = True - except KeyError : - None - - if not trouve : + + elif "GeosCoordinateSystem_h" in self.nc : + # Hrv from msg. + variable = self.nc["GeosCoordinateSystem_h"] + variableX = self.nc["X_h"] + variableY = self.nc["Y_h"] + chaineNavigation = "ImageNavigation_h" + + elif "GeosCoordinateSystem1km" in self.nc : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem1km"] + variableX = self.nc["X1km"] + variableY = self.nc["Y1km"] + chaineNavigation = "ImageNavigation1km" + + elif "GeosCoordinateSystem2km" in self.nc : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem2km"] + variableX = self.nc["X2km"] + variableY = self.nc["Y2km"] + chaineNavigation = "ImageNavigation2km" + + elif "GeosCoordinateSystem" in self.nc : + # Msg in 3 kms. + variable = self.nc["GeosCoordinateSystem"] + variableX = self.nc["X"] + variableY = self.nc["Y"] + chaineNavigation = "ImageNavigation" + + else : print("Variables GeosCoordinateSystemXX not founded.") exit(1) @@ -330,8 +276,8 @@ def res(self) : # print("geotransform = ", geotransform) # geotransform = -5570254, 3000.40604, 0, 5570254, 0, -3000.40604 morceaux = geotransform.split(", ") - self.resolution = float(morceaux[1]) # 3000.40604 - # print("resolution = ", self.resolution) + self.resolution = float(morceaux[1]) + # print("resolution = ", self.resolution) # 3000.40604 self.X = variableX[:] self.nbpix = self.X.shape[0] @@ -345,35 +291,16 @@ def res(self) : self.loff = float(variable.attrs["LOFF"]) # res() - - def get_end_time(self) : - """Get the end time. Global attribute of the netcdf.""" - attr = self.nc.attrs["time_coverage_end"] - # YYYY-MM-DDTHH:MM:SSZNNN - # In some versions milliseconds are present, sometimes not. - # For the goesr : 2024-06-28T10:00:21.1Z - longueurchaine = len(attr) - if longueurchaine == 22 : - # Goesr. - stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") - elif longueurchaine == 20 : - # Mtg. - stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") - else : - # Msg, hima. - stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ%f") - return stacq - # get_end_time() - - def get_start_time(self) : - """Get the start time. Global attribute of the netcdf.""" - attr = self.nc.attrs["time_coverage_start"] - - # YYYY-MM-DDTHH:MM:SSZ + def get_endOrStartTime(self, AttributeName) : + """Get the end or the start time. Global attribute of the netcdf. + AttributName : "time_coverage_start", "time_coverage_end" + """ + attr = self.nc.attrs[AttributeName] + # YYYY-MM-DDTHH:MM:SSZNNN or YYYY-MM-DDTHH:MM:SSZ # In some versions milliseconds are present, sometimes not. longueurchaine = len(attr) if longueurchaine == 22 : - # Goesr. + # Goesr : 2024-06-28T10:00:21.1Z stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") elif longueurchaine == 20 : # Mtg. @@ -382,28 +309,24 @@ def get_start_time(self) : # Msg, hima. stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ%f") return stacq - # get_start_time() + # get_endOrStartTime() @property def start_time(self) : return(self.metadata["start_time"]) + @property def end_time(self) : return(self.metadata["end_time"]) - @property - def alt(self): + def alt(self) : """Get the altitude.""" variable = self.nc["satellite"] - altitude = variable.attrs["dst"] - # 36000000 meters. - altitude += 6378169. # equatorial radius of the earth. + altitude = variable.attrs["dst"] # 36000000 meters. + altitude += 6378169. # equatorial radius of the earth. return altitude - - # ----------------------------------------------------- # - # ----------------------------------------------------- # def preparer_metadata(self, variable) : """Get the metadata.""" mda = {} @@ -413,21 +336,21 @@ def preparer_metadata(self, variable) : mda.update({name: attributs.get(name)}) mda.update({ - "start_time": self.start_time, - "end_time": self.end_time, - "platform_name": self.plateforme, - "sensor": self.sensor, - "zone": self.zone, - "projection_altitude": self.alt, - "cfac": self.cfac, - "lfac": self.lfac, - "coff": self.coff, - "loff": self.loff, - "resolution": self.resolution, - "satellite_actual_longitude": self.longitudeReelle, - "projection_longitude": self.longitudedeprojection, - "projection_type": self.projection - }) + "start_time": self.start_time, + "end_time": self.end_time, + "platform_name": self.plateforme, + "sensor": self.sensor, + "zone": self.zone, + "projection_altitude": self.alt, + "cfac": self.cfac, + "lfac": self.lfac, + "coff": self.coff, + "loff": self.loff, + "resolution": self.resolution, + "satellite_actual_longitude": self.longitudeReelle, + "projection_longitude": self.longitudedeprojection, + "projection_type": self.projection + }) # Placer ici les paramètres orbitaux, une seule fois ? mda.update(self.orbital_param()) @@ -435,143 +358,117 @@ def preparer_metadata(self, variable) : return mda # preparer_metadata(). - def _get_dsname(self, ds_id) : - """Return the correct dataset name based on requested band.""" - # ds_id = DataID(name='vis_08', wavelength=WavelengthRange(min=0.815, central=0.865, max=0.915, unit='µm'), - # resolution=2000, calibration=, modifiers=()) - - # For mtg : - if ds_id["name"] == 'vis_04' : # Name in satpy. - ds_get_name = "VIS004" # Name in icare/meteofrance netcdf. - elif ds_id["name"] == 'vis_05' : - ds_get_name = "VIS005" - elif ds_id["name"] == 'vis_06' : - ds_get_name = "VIS006" - elif ds_id["name"] == 'vis_08' : - ds_get_name = "VIS008" - elif ds_id["name"] == 'vis_09' : - ds_get_name = "VIS009" - elif ds_id["name"] == 'nir_13' : - ds_get_name = "IR_013" - elif ds_id["name"] == 'nir_16' : - ds_get_name = "IR_016" - elif ds_id["name"] == 'nir_22' : - ds_get_name = "IR_022" - elif ds_id["name"] == 'ir_38' : - ds_get_name = "IR_038" - elif ds_id["name"] == 'wv_63' : - ds_get_name = "WV_063" - elif ds_id["name"] == 'wv_73' : - ds_get_name = "WV_073" - elif ds_id["name"] == 'ir_87' : - ds_get_name = "IR_087" - elif ds_id["name"] == 'ir_97' : - ds_get_name = "IR_097" - elif ds_id["name"] == 'ir_105' : - ds_get_name = "IR_105" - elif ds_id["name"] == 'ir_123' : - ds_get_name = "IR_123" - elif ds_id["name"] == 'ir_133' : - ds_get_name = "IR_133" + """Return the correct dataset name based on requested band. + ds_id = DataID(name='vis_08', + wavelength=WavelengthRange(...), + resolution=2000, calibration=, + modifiers=()) + """ + pdict = {} - # For msg, the satpy and icare channel names are the same. - elif ds_id["name"] == 'IR_039' or ds_id["name"] == 'IR_016' or ds_id["name"] == 'VIS008' \ - or ds_id["name"] == 'IR_087' or ds_id["name"] == 'IR_097' or ds_id["name"] == 'IR_108' \ - or ds_id["name"] == 'IR_120' or ds_id["name"] == 'IR_134' or ds_id["name"] == 'VIS006' \ - or ds_id["name"] == 'WV_062' or ds_id["name"] == 'WV_073' or ds_id["name"] == 'HRV' : - ds_get_name = ds_id["name"] - - # For the goesr : - elif ds_id["name"] == 'C01' : - ds_get_name = "VIS_004" - elif ds_id["name"] == 'C02' : - ds_get_name = "VIS_006" - elif ds_id["name"] == 'C03' : - ds_get_name = "VIS_008" - elif ds_id["name"] == 'C04' : - ds_get_name = "VIS_014" - elif ds_id["name"] == 'C05' : - ds_get_name = "VIS_016" - elif ds_id["name"] == 'C06' : - ds_get_name = "VIS_022" - elif ds_id["name"] == 'C07' : - ds_get_name = "IR_039" - elif ds_id["name"] == 'C08' : - ds_get_name = "IR_062" - elif ds_id["name"] == 'C09' : - ds_get_name = "IR_069" - elif ds_id["name"] == 'C10' : - ds_get_name = "IR_073" - elif ds_id["name"] == 'C11' : - ds_get_name = "IR_085" - elif ds_id["name"] == 'C12' : - ds_get_name = "IR_096" - elif ds_id["name"] == 'C13' : - ds_get_name = "IR_103" - elif ds_id["name"] == 'C14' : - ds_get_name = "IR_114" - elif ds_id["name"] == 'C15' : - ds_get_name = "IR_123" - elif ds_id["name"] == 'C16' : - ds_get_name = "IR_133" - - # For himawari. - elif ds_id["name"] == 'B01' : # Name in satpy. - ds_get_name = "VIS004" # Name in icare/meteofrance netcdf. - elif ds_id["name"] == 'B02' : - ds_get_name = "VIS005" - elif ds_id["name"] == 'B03' : - ds_get_name = "VIS006" - elif ds_id["name"] == 'B04' : - ds_get_name = "VIS008" - elif ds_id["name"] == 'B05' : - ds_get_name = "IR_016" - elif ds_id["name"] == 'B06' : - ds_get_name = "IR_022" - elif ds_id["name"] == 'B07' : - ds_get_name = "IR_038" - elif ds_id["name"] == 'B08' : - ds_get_name = "WV_062" - elif ds_id["name"] == 'B09' : - ds_get_name = "WV_069" - elif ds_id["name"] == 'B10' : - ds_get_name = "WV_073" - elif ds_id["name"] == 'B11' : - ds_get_name = "IR_085" - elif ds_id["name"] == 'B12' : - ds_get_name = "IR_096" - elif ds_id["name"] == 'B13' : - ds_get_name = "IR_104" - elif ds_id["name"] == 'B14' : - ds_get_name = "IR_112" - elif ds_id["name"] == 'B15' : - ds_get_name = "IR_123" - elif ds_id["name"] == 'B16' : - ds_get_name = "IR_132" + # For mtg. + # vis_04 is the name in satpy. + # VIS004 is the icare/meteofrance netcdf name. + pdict["vis_04"] = "VIS004" + pdict["vis_05"] = "VIS005" + pdict["vis_06"] = "VIS006" + pdict["vis_08"] = "VIS008" + pdict["vis_09"] = "VIS009" + pdict["nir_13"] = "IR_013" + pdict["nir_16"] = "IR_016" + pdict["nir_22"] = "IR_022" + pdict["ir_38"] = "IR_038" + pdict["wv_63"] = "WV_063" + pdict["wv_73"] = "WV_073" + pdict["ir_87"] = "IR_087" + pdict["ir_97"] = "IR_097" + pdict["ir_105"] = "IR_105" + pdict["ir_123"] = "IR_123" + pdict["ir_133"] = "IR_133" + # For msg, the satpy and icare channel names are the same. + pdict["VIS006"] = "VIS006" + pdict["VIS008"] = "VIS008" + pdict["HRV"] = "HRV" + pdict["IR_016"] = "IR_016" + pdict["IR_039"] = "IR_039" + pdict["WV_062"] = "WV_062" + pdict["WV_073"] = "WV_073" + pdict["IR_087"] = "IR_087" + pdict["IR_097"] = "IR_097" + pdict["IR_108"] = "IR_108" + pdict["IR_120"] = "IR_120" + pdict["IR_134"] = "IR_134" + + # For the goesr satellites : + pdict["C01"] = "VIS_004" + pdict["C02"] = "VIS_006" + pdict["C03"] = "VIS_008" + pdict["C04"] = "VIS_014" + pdict["C05"] = "VIS_016" + pdict["C06"] = "VIS_022" + pdict["C07"] = "IR_039" + pdict["C08"] = "IR_062" + pdict["C09"] = "IR_069" + pdict["C10"] = "IR_073" + pdict["C11"] = "IR_085" + pdict["C12"] = "IR_096" + pdict["C13"] = "IR_103" + pdict["C14"] = "IR_114" + pdict["C15"] = "IR_123" + pdict["C16"] = "IR_133" + + # For himawari. + # BO1 : name in satpy. VIS004 : name in icare/meteofrance netcdf. + pdict["B01"] = "VIS004" + pdict["B02"] = "VIS005" + pdict["B03"] = "VIS006" + pdict["B04"] = "VIS008" + pdict["B05"] = "IR_016" + pdict["B06"] = "IR_022" + pdict["B07"] = "IR_038" + pdict["B08"] = "WV_062" + pdict["B09"] = "WV_069" + pdict["B10"] = "WV_073" + pdict["B11"] = "IR_085" + pdict["B12"] = "IR_096" + pdict["B13"] = "IR_104" + pdict["B14"] = "IR_112" + pdict["B15"] = "IR_123" + pdict["B16"] = "IR_132" + + satpyName = ds_id["name"] + if satpyName in pdict : + icareName = pdict[satpyName] else : - print("Soft not adaptated for this channel : ds_id = ", ds_id["name"]) + print( + "Soft not adaptated for this channel : ds_id = ", + satpyName) exit(1) - return ds_get_name + return icareName # _get_dsname() - def initialisation_dataset(self) : - listeToutesVariables = {"VIS004", "VIS_004", "VIS005", "VIS_005", "VIS006", "VIS_006", "VIS006", \ - "VIS_008", "VIS008", "VIS_009", "VIS009", "IR_013", "IR_016", "IR_022", "IR_038", "IR_039", \ - "IR_062", "WV_062", "WV_063", "IR_069", "WV_069", "IR_073", "WV_073", "IR_085", "IR_087", \ - "IR_096", "IR_097", "IR_103", "IR_104", "IR_105", "IR_108", "IR_112", "IR_120", "IR_123", \ - "IR_132", "IR_133", "IR_134", "HRV"} + listeToutesVariables = { + "VIS004", "VIS_004", "VIS005", "VIS_005", "VIS006", + "VIS_006", "VIS006", "VIS_008", "VIS008", "VIS_009", + "VIS009", "IR_013", "IR_016", "IR_022", "IR_038", + "IR_039", "IR_062", "WV_062", "WV_063", "IR_069", + "WV_069", "IR_073", "WV_073", "IR_085", "IR_087", + "IR_096", "IR_097", "IR_103", "IR_104", "IR_105", + "IR_108", "IR_112", "IR_120", "IR_123", "IR_132", + "IR_133", "IR_134", "HRV" + } self.mda = {} self.scale_factor = {} self.offset = {} self.alpha = {} self.beta = {} - self.bandfactor = {} self.nuc = {} + self.bandfactor = {} self.variableComptes = {} # Loop over the all possible ds_get_name of every satellites : @@ -588,6 +485,7 @@ def initialisation_dataset(self) : # Brightness temperature. self.alpha[ds_get_name] = attributs["alpha"] self.beta[ds_get_name] = attributs["beta"] + self.nuc[ds_get_name] = attributs["nuc"] nomRetourComptes = "Temp_to_Native_count_" + ds_get_name @@ -598,99 +496,132 @@ def initialisation_dataset(self) : nomRetourComptes = "Albedo_to_Native_count_" + ds_get_name else : - print("Nuc or bandfactor not founded int the attributs of ", ds_get_name) + print( + "Nuc or bandfactor not founded int the attributs of ", + ds_get_name) exit(1) self.variableComptes[ds_get_name] = self.nc[nomRetourComptes] - # (65536). Correspondence from 0 to 65535 towards the original spatial agency counts. + # (65536). Correspondence from 0 to 65535 towards + # the original spatial agency counts. self.mda[ds_get_name] = self.preparer_metadata(variable) # initialisation_dataset() - - def get_dataset(self, ds_id, ds_info) : - """Get the dataset. - ds_id["calibration"] = key["calibration"] - = ["brightness_temperature", "reflectance", "radiance", "counts"] - """ - ds_get_name = self._get_dsname(ds_id) # "IR_096" + def comebacktoNativeData(self, ds_get_name) : + """ Come back to the original counts of the hrit. + ds_get_name : meteo france name of a channel : IR_108. """ variable = self.nc[ds_get_name] + # Variable is between -9000 to 4000 (temperature) + # or between 0 to 10000 (albedo). - scale_factor = self.scale_factor[ds_get_name] - offset = self.offset[ds_get_name] + offset = self.offset[ds_get_name] # 0 or 273.15 + variable += 32768 # 0 to 65535 - # print(ds_get_name, ds_id["calibration"], " min = ", np.min(variable[:].values), - # " max = ", np.max(variable[:].values), "dims = ", variable.dims) - # WV_062 calibration.brightness_temperature, from -9000 to 4000 - # VIS006 calibration.reflectance, from 0 to 10000 + if offset == 0. : + # Albedo. + nom = "Albedo_to_Native_count_" + ds_get_name + else : + nom = "Temp_to_Native_count_" + ds_get_name + """ Temp_to_Native_count_IR_062 """ - mda = self.mda[ds_get_name] - mda.update(ds_info) + variableComptes = self.nc[nom] # (65536). + # Correspondence from 0 to 65535 towards the original spatial agency counts. - calibration = ds_id["calibration"] - if calibration == "counts" : - # Come back to the original counts of the hrit... + arrayTableConversion = xr.DataArray.to_numpy(variableComptes) - # Variable is between -9000 to 4000 (temperature) - # or between 0 to 10000 (albedo). + tableau = arrayTableConversion[variable[:]] + """ Come back to the original counts of the hrit. + tableau : 0 to 4095 if native datas coded with 12 bits. """ - variable += 32768 # 0 to 65535 + variable[:] = tableau + return(variable) + # comebacktoNativeData(self, ds_get_name) - if offset == 0. : - # Albedo. - nom = "Albedo_to_Native_count_" + ds_get_name - else : - nom = "Temp_to_Native_count_" + ds_get_name + def comebacktoRadiance(self, ds_get_name) : + # Come back to the radiance. - variableComptes = self.variableComptes[ds_get_name] - # (65536). Correspondence from 0 to 65535 towards the original spatial agency counts. + variable = self.nc[ds_get_name] + # Variable is between -9000 to 4000 (temperature) + # or between 0 to 10000 (albedo). + + scale_factor = self.scale_factor[ds_get_name] # 0.01 + offset = self.offset[ds_get_name] # 0 or 273.15 + if offset == 0. : + # Visible channel. + bandfactor = self.bandfactor[ds_get_name] + + # Variable is an albedo from 0 to 10000. + variable = variable * scale_factor / 100. * bandfactor + # => variable is a reflectance between 0 and 1. + # radiance in mWm-2sr-1(cm-1)-1 + else : + # Brightness temperature. + nuc = self.nuc[ds_get_name] + alpha = self.alpha[ds_get_name] + beta = self.beta[ds_get_name] + + variable = variable * scale_factor + offset + # variable is in Kelvin. + variable = variable * alpha + beta + + c1 = 1.1910427e-5 # Planck + c2 = 1.4387752 # Planck + resul1 = c1 * np.power(nuc, 3.) + resul2 = c2 * nuc + val2 = np.exp(resul2 / variable) - 1. + variable = resul1 / val2 + # Radiance in mWm-2sr-1(cm-1)-1 + return(variable) + # comebacktoRadiance(self, ds_get_name) + + def get_dataset(self, ds_id, ds_info) : + """Get the dataset. + ds_id["calibration"] = key["calibration"] = + ["brightness_temperature", "reflectance", "radiance", "counts"] + """ + ds_get_name = self._get_dsname(ds_id) + # "IR_096" + + mda = self.mda[ds_get_name] + mda.update(ds_info) + + calibration = ds_id["calibration"] + + if calibration == "counts" : # Come back to the original counts of the hrit... - sortie = variablesComptes[AlbedoBT] + variable = self.comebacktoNativeData(ds_get_name) elif calibration == "radiance" : # Come back to the radiance. - - if offset == 0. : - # Visible channel. - bandfactor = self.bandfactor[ds_get_name] - - # Variable is an albedo from 0 to 10000. - variable = variable * scale_factor/ 100. * bandfactor - # => variable is a reflectance between 0 and 1. - # radiance in mWm-2sr-1(cm-1)-1 - else : - # Brightness temperature. - nuc = self.nuc[ds_get_name] - alpha = self.alpha[ds_get_name] - beta = self.beta[ds_get_name] - - variable = variable * scale_factor + offset - # variable is in Kelvin. - variable = variable * alpha + beta - - c1 = 1.1910427e-5 # Planck - c2 = 1.4387752 # Planck - resul1 = c1 * np.power(nuc, 3.) - resul2 = c2 * nuc - val2 = np.exp(resul2 / variable) - 1. - variable = resul1 / val2 - # Radiance in mWm-2sr-1(cm-1)-1 + variable = self.comebacktoRadiance(ds_get_name) elif calibration == "brightness_temperature" : + variable = self.nc[ds_get_name] + # WV_062 calibration.brightness_temperature, from -9000 to 4000 + scale_factor = self.scale_factor[ds_get_name] + offset = self.offset[ds_get_name] if offset != 273.15 : print(ds_get_name, ". offset = ", offset) - print("Soft not intended for a reflectance with a wave length more than 3 microns.") + message = "Soft not intended for a reflectance " + message += "with a wave length more than 3 microns." + print(message) exit(1) variable = variable * scale_factor + offset # variable is in Kelvin. elif calibration == "reflectance" : + variable = self.nc[ds_get_name] + # VIS006 calibration.reflectance, from 0 to 10000 + scale_factor = self.scale_factor[ds_get_name] + offset = self.offset[ds_get_name] if offset != 0. : print(ds_get_name, ". offset = ", offset) - message = "Soft not intended for a brightness temperature " + message = "Soft not intended " + message += "for a brightness temperature " message += "with a wave length less than 3 microns." print(message) exit(1) @@ -707,7 +638,6 @@ def get_dataset(self, ds_id, ds_info) : return variable # get_dataset() - def orbital_param(self) : orb_param_dict = { "orbital_parameters": { @@ -720,10 +650,72 @@ def orbital_param(self) : "projection_longitude": self.longitudedeprojection, "projection_latitude": 0., "projection_altitude": self.alt, - }} + } + } return orb_param_dict + def resolutionSeviri(self, pdict) : + if self.nbpix == 3712 : + pdict["a_desc"] = "MSG/SEVIRI low resolution channel area" + pdict["p_id"] = "msg_lowres" + elif self.nbpix == 11136 : + pdict["a_desc"] = "MSG/SEVIRI HRV channel area" + pdict["p_id"] = "msg_hires" + else : + print("ERROR : not expected resolution for msg : ", self.nbpix) + exit(1) + return(pdict) + + def resolutionFci(self, pdict) : + if self.nbpix == 5568 : + pdict["a_desc"] = "MTG 2km channel area" + pdict["p_id"] = "mtg_lowres" + elif self.nbpix == 11136 : + pdict["a_desc"] = "MTG 1km channel area" + pdict["p_id"] = "mtg_midres" + elif self.nbpix == 22272 : + pdict["a_desc"] = "MTG 500m channel area" + pdict["p_id"] = "mtg_hires" + else : + print("ERROR : not expected resolution for mtg : ", self.nbpix) + exit(1) + return(pdict) + + def resolutionAhi(self, pdict) : + if self.nbpix == 5500 : + pdict["a_desc"] = "HIMA 2km channel area" + pdict["p_id"] = "hima_lowres" + elif self.nbpix == 11000 : + pdict["a_desc"] = "HIMA 1km channel area" + pdict["p_id"] = "hima_midres" + elif self.nbpix == 22000 : + pdict["a_desc"] = "HIMA 500m channel area" + pdict["p_id"] = "hima_hires" + else : + print( + "ERROR : not expected resolution for hima : ", + self.nbpix) + exit(1) + return(pdict) + + def resolutionAbi(self, pdict) : + if self.nbpix == 5424 : + pdict["a_desc"] = "GOESR 2km channel area" + pdict["p_id"] = "goesr_lowres" + elif self.nbpix == 10848 : + pdict["a_desc"] = "GOESR 1km channel area" + pdict["p_id"] = "goesr_midres" + elif self.nbpix == 21696 : + pdict["a_desc"] = "GOESR 500m channel area" + pdict["p_id"] = "goesr_hires" + else : + print( + "ERROR : not expected resolution for goesr : ", + self.nbpix) + exit(1) + return(pdict) + def get_area_def(self, ds_id) : """Get the area def.""" @@ -735,7 +727,8 @@ def get_area_def(self, ds_id) : pdict["a"] = 6378169 pdict["b"] = 6356583.8 - pdict["h"] = self.alt - pdict["a"] # 36000000 mètres. + pdict["h"] = self.alt - pdict["a"] + # 36000000 mètres. pdict["ssp_lon"] = self.longitudedeprojection pdict["ncols"] = self.nblig pdict["nlines"] = self.nbpix @@ -749,67 +742,27 @@ def get_area_def(self, ds_id) : # msg. pdict["scandir"] = "N2S" pdict["a_name"] = "geosmsg" - if self.nbpix == 3712 : - pdict["a_desc"] = "MSG/SEVIRI low resolution channel area" - pdict["p_id"] = "msg_lowres" - elif self.nbpix == 11136 : - pdict["a_desc"] = "MSG/SEVIRI HRV channel area" - pdict["p_id"] = "msg_hires" - else : - print("ERROR : not expected resolution for msg : ", self.nbpix) - exit(1) + pdict = self.resolutionSeviri(pdict) elif self.sensor == "fci" : # mtg. pdict["scandir"] = "N2S" pdict["a_name"] = "geosmtg" - if self.nbpix == 5568 : - pdict["a_desc"] = "MTG 2km channel area" - pdict["p_id"] = "mtg_lowres" - elif self.nbpix == 11136 : - pdict["a_desc"] = "MTG 1km channel area" - pdict["p_id"] = "mtg_midres" - elif self.nbpix == 22272 : - pdict["a_desc"] = "MTG 500m channel area" - pdict["p_id"] = "mtg_hires" - else : - print("ERROR : not expected resolution for mtg : ", self.nbpix) - exit(1) + pdict = self.resolutionFci(pdict) elif self.sensor == "ahi" : # Himawari. pdict["scandir"] = "N2S" pdict["a_name"] = "geoshima" - if self.nbpix == 5500 : - pdict["a_desc"] = "HIMA 2km channel area" - pdict["p_id"] = "hima_lowres" - elif self.nbpix == 11000 : - pdict["a_desc"] = "HIMA 1km channel area" - pdict["p_id"] = "hima_midres" - elif self.nbpix == 22000 : - pdict["a_desc"] = "HIMA 500m channel area" - pdict["p_id"] = "hima_hires" - else : - print("ERROR : not expected resolution for hima : ", self.nbpix) - exit(1) + pdict = self.resolutionAhi(pdict) elif self.sensor == "abi" : # Goesr. pdict["scandir"] = "N2S" pdict["a_name"] = "geosgoesr" pdict["sweep"] = "x" - if self.nbpix == 5424 : - pdict["a_desc"] = "GOESR 2km channel area" - pdict["p_id"] = "goesr_lowres" - elif self.nbpix == 10848 : - pdict["a_desc"] = "GOESR 1km channel area" - pdict["p_id"] = "goesr_midres" - elif self.nbpix == 21696 : - pdict["a_desc"] = "GOESR 500m channel area" - pdict["p_id"] = "goesr_hires" - else : - print("ERROR : not expected resolution for goesr : ", self.nbpix) - exit(1) + pdict = self.resolutionAbi(pdict) + else : print("ERROR : " + self.sensor + " not expected.") exit(1) diff --git a/satpy/tests/reader_tests/test_geos_netcdficare.py b/satpy/tests/reader_tests/test_geos_netcdficare.py new file mode 100644 index 0000000000..cf1f5401ad --- /dev/null +++ b/satpy/tests/reader_tests/test_geos_netcdficare.py @@ -0,0 +1,308 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Copyright (c) 2019 Satpy developers +# +# This file is part of satpy. +# +# satpy is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# satpy is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along with +# satpy. If not, see . + +"""Tests for the Icare MeteoFrance netcdfs reader, + satpy/readers/geos_netcdficare.py. + +SYNTAXE : + pytest + + Before, copy the satpy/tests/reader_tests/test_geos_netcdfcare.py file + into the pytest directory pytest_projec. + +EXECUTION TIME : + 20 seconds. +DATE OF CREATION : + 2024 11th october. +LAST VERSIONS : + +AUTHOR : + Meteo France. +""" + +import os + +import numpy as np + +from satpy.scene import Scene +from satpy import find_files_and_readers + +from datetime import datetime +from netCDF4 import Dataset + + +class TestGeosNetcdfIcareReader() : + # Test of the geos_netcdficare reader. + # This reader has been build for the Icare Meteo France netcdfs. + + def test_geos_netcdficare(self, tmp_path) : + """ A dummy netcdf is built. + A scene self.scn for the convection product for the same date + is built. We check that the scene parameters are the same + as thoses in the dummy netcdf. + This procedure is called by pytest. + """ + + self.init(tmp_path) + + startTime = self.scn.start_time + startTimeString = startTime.strftime('%Y-%m-%dT%H:%M:%S') + # 2024-06-28T10:00:40 + assert startTimeString == self.expectedStartTime + + endTime = self.scn.end_time + endTimeString = endTime.strftime('%Y-%m-%dT%H:%M:%S') + # 2024-06-28T10:12:41 + assert endTimeString == self.expectedEndTime + + sensor = self.scn.sensor_names + for isensor in sensor : + capteur = isensor + assert capteur == self.expectedSensor + + platform = "erreur" + altitude = -1. + longitude = 999. + + for data_arr in self.values : + # values come from the scene. + # print("data_arr.attrs = ", data_arr.attrs) + if "platform_name" in data_arr.attrs : + platform = data_arr.attrs["platform_name"] + if "orbital_parameters" in data_arr.attrs : + subAttr = data_arr.attrs["orbital_parameters"] + if "satellite_actual_altitude" in subAttr : + altitude = subAttr["satellite_actual_altitude"] + if "satellite_actual_longitude" in data_arr.attrs : + longitude = data_arr.attrs["satellite_actual_longitude"] + + longitude = float(int(longitude * 10.)) / 10. + assert platform == self.expectedPlatform + assert longitude == self.expectedLongitude + assert altitude == self.expectedAltitude + + xr = self.scn.to_xarray_dataset() + # print("xr = ", xr) + matrice = xr["convection"] + nblin = matrice.shape[1] + nbpix = matrice.shape[2] + assert nblin == self.expectedNblin + assert nbpix == self.expectedNbpix + + cfac = xr.attrs["cfac"] + assert cfac == self.expectedCfac + lfac = xr.attrs["lfac"] + assert lfac == self.expectedLfac + coff = xr.attrs["coff"] + assert coff == self.expectedCoff + loff = xr.attrs["loff"] + assert loff == self.expectedLoff + + satpyId = xr.attrs["_satpy_id"] + # DataID(name='convection', resolution=3000.403165817) + # Cf satpy/dataset/dataid.py. + + resolution = satpyId.get("resolution") + resolution = float(int(resolution * 10.)) / 10. + assert resolution == self.expectedResolution + + # A picture of convection composite will be displayed. + self.scn.show("convection") + print("The picture should be pink.") + # test_geos_netcdficare(self, tmp_path) + + def init(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_geos_netcdficare(). + """ + self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" + self.filepath = tmp_path + + self.buildNetcdf(self.netcdfName) + + # We will check that the parameters written in the dummy netcdf can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:12:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "Meteosat-10" + self.expectedSensor = "seviri" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 1.3642337E7 + self.expectedLfac = 1.3642337E7 + self.expectedCoff = 1857.0 + self.expectedLoff = 1857.0 + self.expectedResolution = 3000.4 + self.expectedNblin = 3712 + self.expectedNbpix = 3712 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for in the filepath directory. + yaml_file = 'msg_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + print("Found myfiles = ", myfiles) + # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + + print(self.scn.available_dataset_names()) + # ['IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', + # 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073'] + + print(self.scn.available_composite_names()) + """ Static compositor {'_satpy_id': DataID(name='_night_background'), + 'standard_name': 'night_background', + 'prerequisites': [], + 'optional_prerequisites': []} + Static compositor {'_satpy_id': DataID(name='_night_background_hires'), + 'standard_name': 'night_background', 'prerequisites': [], + 'optional_prerequisites': []} + ['airmass', 'ash', 'cloud_phase_distinction', 'cloud_phase_distinction_raw', + 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection'... + """ + + self.scn.load(['convection']) + self.values = self.scn.values() + # init() + + def buildNetcdf(self, ncName) : + """ + ncName : tmp_path / Mmultic3kmNC4_msg03_202406281000.nc + A dummy icare Meteo France netcdf is built here. + Called by init(). + """ + if os.path.exists(ncName) : + os.remove(ncName) + ncfileOut = Dataset( + ncName, mode="w", clobber=True, + diskless=False, persist=False, format='NETCDF4') + + ncfileOut.createDimension(u'ny', 3712) + ncfileOut.createDimension(u'nx', 3712) + ncfileOut.createDimension(u'numerical_count', 65536) + ncfileOut.setncattr("time_coverage_start", "2024-06-28T10:00:09Z383") + ncfileOut.setncattr("time_coverage_end", "2024-06-28T10:12:41Z365") + ncfileOut.setncattr("Area_of_acquisition", "globe") + + fill_value = -32768 + var = ncfileOut.createVariable( + "satellite", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + + var.setncattr("id", "msg03") + var.setncattr("dst", 35786691.) + var.setncattr("lon", float(0.1)) + + var = ncfileOut.createVariable( + "geos", "c", zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None) + var.setncattr("longitude_of_projection_origin", 0.) + + var = ncfileOut.createVariable( + "GeosCoordinateSystem", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + var.setncattr( + "GeoTransform", + "-5570254, 3000.40604, 0, 5570254, 0, -3000.40604") + + var = ncfileOut.createVariable( + "ImageNavigation", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + var.setncattr("CFAC", 1.3642337E7) + var.setncattr("LFAC", 1.3642337E7) + var.setncattr("COFF", 1857.0) + var.setncattr("LOFF", 1857.0) + + var = ncfileOut.createVariable( + "X", 'float32', u'nx', zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + x0 = -5570254. + dx = 3000.40604 + var[:] = np.array(([(x0 + dx * i) for i in range(3712)])) + + var = ncfileOut.createVariable( + "Y", 'float32', u'ny', zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + y0 = 5570254. + dy = -3000.40604 + var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) + + for channel in {"VIS006", "VIS008", "IR_016"} : + var = ncfileOut.createVariable( + channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', + least_significant_digit=None, fill_value=fill_value) + var[:] = np.array(([[i * 2 for i in range(3712)] for j in range(3712)])) + # Hundredths of albedo between 0 and 10000. + var.setncattr("scale_factor", 0.01) + var.setncattr("add_offset", 0.) + var.setncattr("bandfactor", 20.76) + var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + + var = ncfileOut.createVariable( + "Albedo_to_Native_count_" + channel, 'short', + 'numerical_count', zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None, + fill_value=-9999) + var[:] = np.array(([-9999 for i in range(65536)])) + # In order to come back to the native datas on 10, 12 or 16 bits. + + for channel in { + "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", + "IR_108", "IR_120", "IR_134"} : + var = ncfileOut.createVariable( + channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', + least_significant_digit=None, fill_value=fill_value) + var[:] = np.array( + ([[-9000 + j * 4 for i in range(3712)] for j in range(3712)])) + # Hundredths of celcius degrees. + var.setncattr("scale_factor", 0.01) + var.setncattr("add_offset", 273.15) + var.setncattr("nuc", 1600.548) + var.setncattr("alpha", 0.9963) + var.setncattr("beta", 2.185) + var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + + var = ncfileOut.createVariable( + "Temp_to_Native_count_" + channel, 'short', + 'numerical_count', zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None, + fill_value=-9999) + var[:] = np.array(([-9999 for i in range(65536)])) + # In order to come back to the native datas on 10, 12 or 16 bits. + ncfileOut.close + # buildNetcdf() + + # class TestGeosNetcdfIcareReader. From b07d8a8277f8d5f3b62d8489173c1db3aac510ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?laurent=20p=C3=A9rier?= Date: Tue, 15 Oct 2024 11:39:55 +0200 Subject: [PATCH 5/9] Answer to codescene request. --- satpy/readers/geos_netcdficare.py | 160 +++++++++--------- .../reader_tests/test_geos_netcdficare.py | 20 ++- 2 files changed, 90 insertions(+), 90 deletions(-) diff --git a/satpy/readers/geos_netcdficare.py b/satpy/readers/geos_netcdficare.py index d4cd1cff12..154aa870f0 100644 --- a/satpy/readers/geos_netcdficare.py +++ b/satpy/readers/geos_netcdficare.py @@ -155,19 +155,24 @@ def sensor_name(self) : The sensor and platform names are stored together, eg: MSG1/SEVIRI """ variable = self.nc["satellite"] - self.plateforme = variable.attrs["id"] - - if "goes" in self.plateforme : - self.sensor = "abi" - elif ("msg" in self.plateforme) or ("MSG" in self.plateforme) : - self.sensor = "seviri" - elif ("mtg" in self.plateforme) or ("MTG" in self.plateforme) : - self.sensor = "fci" - elif "hima" in self.plateforme : - self.sensor = "ahi" + self.plateforme = variable.attrs["id"] # msg1, msg01, MSG1... + platform = self.plateforme[:3] # msg, MSG + platform = platform.lower() # msg + + pdict = {} + pdict["msg"] = "seviri" + pdict["mtg"] = "fci" + pdict["goe"] = "abi" + pdict["him"] = "ahi" + + if platform in pdict : + self.sensor = pdict[platform] else : - raise NameError("Unsupported satellite platform : " + self.plateforme) - + print( + "Unsupported satellite platform : " + + self.plateforme) + exit(1) + # Icare and météo france use non-standard platform names. # Change is needed for pyspectral : # pyspectral/rsr_seviri_Meteosat-10.h5 in the call @@ -358,13 +363,7 @@ def preparer_metadata(self, variable) : return mda # preparer_metadata(). - def _get_dsname(self, ds_id) : - """Return the correct dataset name based on requested band. - ds_id = DataID(name='vis_08', - wavelength=WavelengthRange(...), - resolution=2000, calibration=, - modifiers=()) - """ + def buildChannelCorrespondanceName(self) : pdict = {} # For mtg. @@ -419,7 +418,7 @@ def _get_dsname(self, ds_id) : pdict["C15"] = "IR_123" pdict["C16"] = "IR_133" - # For himawari. + # For himawari. # BO1 : name in satpy. VIS004 : name in icare/meteofrance netcdf. pdict["B01"] = "VIS004" pdict["B02"] = "VIS005" @@ -437,13 +436,24 @@ def _get_dsname(self, ds_id) : pdict["B14"] = "IR_112" pdict["B15"] = "IR_123" pdict["B16"] = "IR_132" + return pdict + # buildChannelCorrespondanceName() + + def _get_dsname(self, ds_id) : + """Return the correct dataset name based on requested band. + ds_id = DataID(name='vis_08', + wavelength=WavelengthRange(...), + resolution=2000, calibration=, + modifiers=()) + """ + pdict = self.buildChannelCorrespondanceName() satpyName = ds_id["name"] if satpyName in pdict : icareName = pdict[satpyName] else : print( - "Soft not adaptated for this channel : ds_id = ", + "Soft not adaptated for this channel : ds_id = ", satpyName) exit(1) @@ -484,8 +494,8 @@ def initialisation_dataset(self) : if "nuc" in attributs : # Brightness temperature. self.alpha[ds_get_name] = attributs["alpha"] - self.beta[ds_get_name] = attributs["beta"] - self.nuc[ds_get_name] = attributs["nuc"] + self.beta[ds_get_name] = attributs["beta"] + self.nuc[ds_get_name] = attributs["nuc"] nomRetourComptes = "Temp_to_Native_count_" + ds_get_name @@ -655,64 +665,15 @@ def orbital_param(self) : return orb_param_dict - def resolutionSeviri(self, pdict) : - if self.nbpix == 3712 : - pdict["a_desc"] = "MSG/SEVIRI low resolution channel area" - pdict["p_id"] = "msg_lowres" - elif self.nbpix == 11136 : - pdict["a_desc"] = "MSG/SEVIRI HRV channel area" - pdict["p_id"] = "msg_hires" - else : - print("ERROR : not expected resolution for msg : ", self.nbpix) - exit(1) - return(pdict) - - def resolutionFci(self, pdict) : - if self.nbpix == 5568 : - pdict["a_desc"] = "MTG 2km channel area" - pdict["p_id"] = "mtg_lowres" - elif self.nbpix == 11136 : - pdict["a_desc"] = "MTG 1km channel area" - pdict["p_id"] = "mtg_midres" - elif self.nbpix == 22272 : - pdict["a_desc"] = "MTG 500m channel area" - pdict["p_id"] = "mtg_hires" - else : - print("ERROR : not expected resolution for mtg : ", self.nbpix) - exit(1) - return(pdict) - - def resolutionAhi(self, pdict) : - if self.nbpix == 5500 : - pdict["a_desc"] = "HIMA 2km channel area" - pdict["p_id"] = "hima_lowres" - elif self.nbpix == 11000 : - pdict["a_desc"] = "HIMA 1km channel area" - pdict["p_id"] = "hima_midres" - elif self.nbpix == 22000 : - pdict["a_desc"] = "HIMA 500m channel area" - pdict["p_id"] = "hima_hires" - else : - print( - "ERROR : not expected resolution for hima : ", - self.nbpix) - exit(1) - return(pdict) - - def resolutionAbi(self, pdict) : - if self.nbpix == 5424 : - pdict["a_desc"] = "GOESR 2km channel area" - pdict["p_id"] = "goesr_lowres" - elif self.nbpix == 10848 : - pdict["a_desc"] = "GOESR 1km channel area" - pdict["p_id"] = "goesr_midres" - elif self.nbpix == 21696 : - pdict["a_desc"] = "GOESR 500m channel area" - pdict["p_id"] = "goesr_hires" + def channelType(self, pdict, pdictResoAdesc, pdictResoPid, satellite) : + strNbpix = str(self.nbpix) + if strNbpix in pdictResoAdesc : + pdict["a_desc"] = pdictResoAdesc[strNbpix] + pdict["p_id"] = pdictResoPid[strNbpix] else : print( - "ERROR : not expected resolution for goesr : ", - self.nbpix) + "ERROR : not expected resolution for ", + satellite, " : ", self.nbpix) exit(1) return(pdict) @@ -738,30 +699,63 @@ def get_area_def(self, ds_id) : pdict["scandir"] = "S2N" pdict["a_name"] = "geosmsg" + pdictResoAdesc = {} + pdictResoPid = {} + if self.sensor == "seviri" : # msg. pdict["scandir"] = "N2S" pdict["a_name"] = "geosmsg" - pdict = self.resolutionSeviri(pdict) + + pdictResoAdesc["3712"] = "MSG/SEVIRI low resolution channel area" + pdictResoPid["3712"] = "msg_lowres" + pdictResoAdesc["11136"] = "MSG/SEVIRI HRV channel area" + pdictResoPid["11136"] = "msg_hires" + + pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "msg") elif self.sensor == "fci" : # mtg. pdict["scandir"] = "N2S" pdict["a_name"] = "geosmtg" - pdict = self.resolutionFci(pdict) + + pdictResoAdesc["5568"] = "MTG 2km channel area" + pdictResoPid["5568"] = "mtg_lowres" + pdictResoAdesc["11136"] = "MTG 1km channel area" + pdictResoPid["11136"] = "mtg_midres" + pdictResoAdesc["22272"] = "MTG 500m channel area" + pdictResoPid["22272"] = "mtg_hires" + + pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "mtg") elif self.sensor == "ahi" : # Himawari. pdict["scandir"] = "N2S" pdict["a_name"] = "geoshima" - pdict = self.resolutionAhi(pdict) + + pdictResoAdesc["5500"] = "HIMA 2km channel area" + pdictResoPid["5500"] = "hima_lowres" + pdictResoAdesc["11000"] = "HIMA 1km channel area" + pdictResoPid["11000"] = "hima_midres" + pdictResoAdesc["22000"] = "HIMA 500m channel area" + pdictResoPid["22000"] = "hima_hires" + + pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "hima") elif self.sensor == "abi" : # Goesr. pdict["scandir"] = "N2S" pdict["a_name"] = "geosgoesr" pdict["sweep"] = "x" - pdict = self.resolutionAbi(pdict) + + pdictResoAdesc["5424"] = "GOESR 2km channel area" + pdictResoPid["5424"] = "goesr_lowres" + pdictResoAdesc["10848"] = "GOESR 1km channel area" + pdictResoPid["10848"] = "goesr_midres" + pdictResoAdesc["21696"] = "GOESR 500m channel area" + pdictResoPid["21696"] = "goesr_hires" + + pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "goesr") else : print("ERROR : " + self.sensor + " not expected.") diff --git a/satpy/tests/reader_tests/test_geos_netcdficare.py b/satpy/tests/reader_tests/test_geos_netcdficare.py index cf1f5401ad..2cda858b52 100644 --- a/satpy/tests/reader_tests/test_geos_netcdficare.py +++ b/satpy/tests/reader_tests/test_geos_netcdficare.py @@ -50,10 +50,10 @@ class TestGeosNetcdfIcareReader() : # Test of the geos_netcdficare reader. # This reader has been build for the Icare Meteo France netcdfs. - def test_geos_netcdficare(self, tmp_path) : - """ A dummy netcdf is built. + def test_geosNetcdficare(self, tmp_path) : + """ A dummy netcdf is built. A scene self.scn for the convection product for the same date - is built. We check that the scene parameters are the same + is built. We check that the scene parameters are the same as thoses in the dummy netcdf. This procedure is called by pytest. """ @@ -124,7 +124,7 @@ def test_geos_netcdficare(self, tmp_path) : # A picture of convection composite will be displayed. self.scn.show("convection") print("The picture should be pink.") - # test_geos_netcdficare(self, tmp_path) + # test_geosNetcdficare(self, tmp_path) def init(self, tmp_path) : """ @@ -155,7 +155,7 @@ def init(self, tmp_path) : self.expectedNbpix = 3712 # To build a scene at the date 20240628_100000, - # a netcdf corresponding to msg_netcdficare + # a netcdf corresponding to msg_netcdficare # is looked for in the filepath directory. yaml_file = 'msg_netcdficare' myfiles = find_files_and_readers( @@ -254,6 +254,13 @@ def buildNetcdf(self, ncName) : dy = -3000.40604 var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) + self.visibleChannelsCreation(ncfileOut, fill_value) + self.infrarougeChannelsCreation(ncfileOut, fill_value) + + ncfileOut.close + # buildNetcdf() + + def visibleChannelsCreation(self, ncfileOut, fill_value) : for channel in {"VIS006", "VIS008", "IR_016"} : var = ncfileOut.createVariable( channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, @@ -276,6 +283,7 @@ def buildNetcdf(self, ncName) : var[:] = np.array(([-9999 for i in range(65536)])) # In order to come back to the native datas on 10, 12 or 16 bits. + def infrarougeChannelsCreation(self, ncfileOut, fill_value) : for channel in { "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", "IR_108", "IR_120", "IR_134"} : @@ -302,7 +310,5 @@ def buildNetcdf(self, ncName) : fill_value=-9999) var[:] = np.array(([-9999 for i in range(65536)])) # In order to come back to the native datas on 10, 12 or 16 bits. - ncfileOut.close - # buildNetcdf() # class TestGeosNetcdfIcareReader. From e4d53b77ff3d253c94b290baf1cd94396a2c2bec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?laurent=20p=C3=A9rier?= Date: Tue, 15 Oct 2024 12:49:06 +0200 Subject: [PATCH 6/9] =?UTF-8?q?Answer=20to=20codescene=20request=20n=C2=B0?= =?UTF-8?q?3.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- satpy/etc/readers/msg_netcdficare.yaml | 14 ++++++++------ satpy/tests/reader_tests/test_geos_netcdficare.py | 2 +- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/satpy/etc/readers/msg_netcdficare.yaml b/satpy/etc/readers/msg_netcdficare.yaml index 8885625321..8a64c7ed01 100644 --- a/satpy/etc/readers/msg_netcdficare.yaml +++ b/satpy/etc/readers/msg_netcdficare.yaml @@ -29,12 +29,12 @@ file_types: MSG1km : file_reader: !!python/name:satpy.readers.geos_netcdficare.NETCDF_ICARE - file_patterns: [ 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', - 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', - 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', - 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', - 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', - 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc'] + file_patterns: ['Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d%H%M}.nc', + 'Mmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Imultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc', + 'Mrsmultic1kmNC4_{platform_shortname:5s}_{start_time:%Y%m%d_%H%M}.nc'] # platform_shortname:5s : msg01, ..., msg04 datasets: @@ -170,3 +170,5 @@ datasets: standard_name: brightness_temperature units: "K" file_type: MSG3km + + diff --git a/satpy/tests/reader_tests/test_geos_netcdficare.py b/satpy/tests/reader_tests/test_geos_netcdficare.py index 2cda858b52..75e27e6dd7 100644 --- a/satpy/tests/reader_tests/test_geos_netcdficare.py +++ b/satpy/tests/reader_tests/test_geos_netcdficare.py @@ -130,7 +130,7 @@ def init(self, tmp_path) : """ A fake netcdf is built. A scene is built with the reader to be tested, applied to this netcdf. - Called by test_geos_netcdficare(). + Called by test_geosNetcdficare(). """ self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" self.filepath = tmp_path From b515ec322c9f1dc1e34d6ebc314eae1308655a24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?laurent=20p=C3=A9rier?= Date: Wed, 16 Oct 2024 16:42:49 +0200 Subject: [PATCH 7/9] Coding guidelines recommandations. --- satpy/readers/geos_netcdficare.py | 1360 ++++++++--------- .../reader_tests/test_geos_netcdficare.py | 548 +++---- test_geos_netcdficare.py | 308 ---- 3 files changed, 954 insertions(+), 1262 deletions(-) delete mode 100644 test_geos_netcdficare.py diff --git a/satpy/readers/geos_netcdficare.py b/satpy/readers/geos_netcdficare.py index 154aa870f0..39fbef1030 100644 --- a/satpy/readers/geos_netcdficare.py +++ b/satpy/readers/geos_netcdficare.py @@ -24,7 +24,7 @@ The ``geos_netcdficare`` reader reads some geostationnary netcdf build by Meteo France and stored at Icare. -The brightness tempera ture and albedo are calibrated. +The brightness temperature and albedo are calibrated. That has been stored by the ICARE Data and Services Center Data can be accessed via: http://www.icare.univ-lille1.fr @@ -67,23 +67,32 @@ -------- Here is an example how to read the data in satpy: - from satpy import Scene - import glob + from satpy import Scene + import glob - filenames = glob.glob('data/*2019-03-01T12-00-00*.hdf') - scn = Scene(filenames = filenames, reader = 'hima_netcdficare') - scn.load(['true_color']) # scn.load(['VIS006']) + filenames = glob.glob('data/*2019-03-01T12-00-00*.hdf') + scn = Scene(filenames = filenames, reader = 'hima_netcdficare') + scn.load(['true_color']) # scn.load(['VIS006']) - my_area = AreaDefinition( - 'my_area', 'zone', 'my_area', - '+proj=latlong +lon_0=0 +a=6378169 +b=6356583 +h=35785831 +x_0=0 - +y_0=0 +pm=0', - 8500, 4000, - [-180., -80., 180., 80], - nprocs=16) + my_area = AreaDefinition( + 'my_area', 'zone', 'my_area', + '+proj=latlong +lon_0=0 +a=6378169 +b=6356583 +h=35785831 +x_0=0 + +y_0=0 +pm=0', + 8500, 4000, + [-180., -80., 180., 80], + nprocs=16) - natscn = scn.resample(my_area, resampler='nearest') - natscn.save_dataset(composite_name, filename = filename_image_out) + natscn = scn.resample(my_area, resampler='nearest') + natscn.save_dataset(composite_name, filename = filename_image_out) + +EXECUTION TIME : + 50 seconds for a 2 kms goesr airmass rgb disk. +DATE OF CREATION : + 2024 16th october. +LAST VERSIONS : + +AUTHOR : + Meteo France. """ @@ -93,677 +102,664 @@ from satpy.readers._geos_area import get_area_definition, get_area_extent -# from netCDF4 import Dataset -# from xarray import DataArray import xarray as xr from satpy.readers.file_handlers import BaseFileHandler -# from satpy._compat import cached_property -# from satpy.utils import get_dask_chunk_size_in_bytes -import dask -from satpy.readers import open_file_or_filename -# import math +import logging +logger = logging.getLogger('netcdficare') + +# Planck : +C1 = 1.1910427e-5 +C2 = 1.4387752 class NETCDF_ICARE(BaseFileHandler) : - # Cf readers/file_handlers.py. - - def __init__(self, filename, filename_info, filetype_info) : - """Init the file handler.""" - - super().__init__(filename, filename_info, filetype_info) - - # chunk_bytes = self._chunk_bytes_for_resolution() - chunk_bytes = '128 MiB' - # chunk_bytes = '64 MiB' - - # with dask.config.set(**{'array.slicing.split_large_chunks': True}) : - with dask.config.set({"array.chunk-size": chunk_bytes}) : - f_obj = open_file_or_filename(self.filename) - self.nc = xr.open_dataset( - f_obj, decode_cf=True, mask_and_scale=False, - chunks={"xc": "auto", "yc": "auto"}) - - self.metadata = {} - - self.metadata["start_time"] = self.get_endOrStartTime( - "time_coverage_start") - self.metadata["end_time"] = self.get_endOrStartTime( - "time_coverage_end") - - print( - "Reading: {}".format(filename), - " start: {} ".format(self.start_time), - " end: {}".format(self.end_time)) - self._cache = {} - - self.sensor_name() - self.res() - self.longitudeReelle = self.satlon() - self.longitudedeprojection = self.projlon() - - self.zone = self.nc.attrs["Area_of_acquisition"] - # globe, europe. - - # Reading the needed attributes. - self.initialisation_dataset() - # __init__() - - def sensor_name(self) : - """Get the sensor name. - The sensor and platform names are stored together, eg: MSG1/SEVIRI - """ - variable = self.nc["satellite"] - self.plateforme = variable.attrs["id"] # msg1, msg01, MSG1... - platform = self.plateforme[:3] # msg, MSG - platform = platform.lower() # msg - - pdict = {} - pdict["msg"] = "seviri" - pdict["mtg"] = "fci" - pdict["goe"] = "abi" - pdict["him"] = "ahi" - - if platform in pdict : - self.sensor = pdict[platform] - else : - print( - "Unsupported satellite platform : " + - self.plateforme) - exit(1) - - # Icare and météo france use non-standard platform names. - # Change is needed for pyspectral : - # pyspectral/rsr_seviri_Meteosat-10.h5 in the call - # Calculator(platform_name, sensor, name). - pdict = {} - pdict["msg1"] = "Meteosat-08" - pdict["msg01"] = "Meteosat-08" - pdict["MSG1"] = "Meteosat-08" - pdict["msg2"] = "Meteosat-09" - pdict["msg02"] = "Meteosat-09" - pdict["MSG2"] = "Meteosat-09" - pdict["msg3"] = "Meteosat-10" - pdict["msg03"] = "Meteosat-10" - pdict["MSG3"] = "Meteosat-10" - pdict["msg4"] = "Meteosat-11" - pdict["msg04"] = "Meteosat-11" - pdict["MSG4"] = "Meteosat-11" - pdict["mtgi1"] = "Meteosat-12" - pdict["mtg1"] = "Meteosat-12" - pdict["MTG01"] = "Meteosat-12" - pdict["goes16"] = "GOES-16" - pdict["goes17"] = "GOES-17" - pdict["goes18"] = "GOES-18" - pdict["goes19"] = "GOES-19" - pdict["hima08"] = "Himawari-8" - pdict["hima09"] = "Himawari-9" - - if self.plateforme in pdict : - self.plateforme = pdict[self.plateforme] - else : - print( - "Unsupported satellite platform : " + - self.plateforme) - exit(1) - - # sensor_name() - - def satlon(self) : - """Get the satellite longitude.""" - variable = self.nc["satellite"] - longitudeReelle = variable.attrs["lon"] - return longitudeReelle - - def projlon(self): - """Get the projection longitude.""" - variable = self.nc["geos"] - longitudedeprojection = variable.attrs["longitude_of_projection_origin"] - return longitudedeprojection - - @property - def projection(self): - """Get the projection.""" - return "geos" - - def res(self) : - """Get the resolution. - The resolution can be read in the attribute geotransform - of the following variables : - GeosCoordinateSystem500m, GeosCoordinateSystem_h, - GeosCoordinateSystem1km, GeosCoordinateSystem2km, - GeosCoordinateSystem. - cfac, lfac, coff, loff can be read in the variables ImageNavigationxxx. - """ - - if "GeosCoordinateSystem500m" in self.nc : - # Mtg, himawari, goesr. - variable = self.nc["GeosCoordinateSystem500m"] - variableX = self.nc["X500m"] - variableY = self.nc["Y500m"] - chaineNavigation = "ImageNavigation500m" - - elif "GeosCoordinateSystem_h" in self.nc : - # Hrv from msg. - variable = self.nc["GeosCoordinateSystem_h"] - variableX = self.nc["X_h"] - variableY = self.nc["Y_h"] - chaineNavigation = "ImageNavigation_h" - - elif "GeosCoordinateSystem1km" in self.nc : - # Mtg, himawari, goesr. - variable = self.nc["GeosCoordinateSystem1km"] - variableX = self.nc["X1km"] - variableY = self.nc["Y1km"] - chaineNavigation = "ImageNavigation1km" - - elif "GeosCoordinateSystem2km" in self.nc : - # Mtg, himawari, goesr. - variable = self.nc["GeosCoordinateSystem2km"] - variableX = self.nc["X2km"] - variableY = self.nc["Y2km"] - chaineNavigation = "ImageNavigation2km" - - elif "GeosCoordinateSystem" in self.nc : - # Msg in 3 kms. - variable = self.nc["GeosCoordinateSystem"] - variableX = self.nc["X"] - variableY = self.nc["Y"] - chaineNavigation = "ImageNavigation" - - else : - print("Variables GeosCoordinateSystemXX not founded.") - exit(1) - - geotransform = variable.attrs["GeoTransform"] - - # print("geotransform = ", geotransform) - # geotransform = -5570254, 3000.40604, 0, 5570254, 0, -3000.40604 - morceaux = geotransform.split(", ") - self.resolution = float(morceaux[1]) - # print("resolution = ", self.resolution) # 3000.40604 - - self.X = variableX[:] - self.nbpix = self.X.shape[0] - self.Y = variableY[:] - self.nblig = self.Y.shape[0] - - variable = self.nc[chaineNavigation] - self.cfac = float(variable.attrs["CFAC"]) - self.lfac = float(variable.attrs["LFAC"]) - self.coff = float(variable.attrs["COFF"]) - self.loff = float(variable.attrs["LOFF"]) - # res() - - def get_endOrStartTime(self, AttributeName) : - """Get the end or the start time. Global attribute of the netcdf. - AttributName : "time_coverage_start", "time_coverage_end" - """ - attr = self.nc.attrs[AttributeName] - # YYYY-MM-DDTHH:MM:SSZNNN or YYYY-MM-DDTHH:MM:SSZ - # In some versions milliseconds are present, sometimes not. - longueurchaine = len(attr) - if longueurchaine == 22 : - # Goesr : 2024-06-28T10:00:21.1Z - stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") - elif longueurchaine == 20 : - # Mtg. - stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") - else : - # Msg, hima. - stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ%f") - return stacq - # get_endOrStartTime() - - @property - def start_time(self) : - return(self.metadata["start_time"]) - - @property - def end_time(self) : - return(self.metadata["end_time"]) - - @property - def alt(self) : - """Get the altitude.""" - variable = self.nc["satellite"] - altitude = variable.attrs["dst"] # 36000000 meters. - altitude += 6378169. # equatorial radius of the earth. - return altitude - - def preparer_metadata(self, variable) : - """Get the metadata.""" - mda = {} - - attributs = variable.attrs - for name in attributs : - mda.update({name: attributs.get(name)}) - - mda.update({ - "start_time": self.start_time, - "end_time": self.end_time, - "platform_name": self.plateforme, - "sensor": self.sensor, - "zone": self.zone, - "projection_altitude": self.alt, - "cfac": self.cfac, - "lfac": self.lfac, - "coff": self.coff, - "loff": self.loff, - "resolution": self.resolution, - "satellite_actual_longitude": self.longitudeReelle, - "projection_longitude": self.longitudedeprojection, - "projection_type": self.projection - }) - - # Placer ici les paramètres orbitaux, une seule fois ? - mda.update(self.orbital_param()) - - return mda - # preparer_metadata(). - - def buildChannelCorrespondanceName(self) : - pdict = {} - - # For mtg. - # vis_04 is the name in satpy. - # VIS004 is the icare/meteofrance netcdf name. - pdict["vis_04"] = "VIS004" - pdict["vis_05"] = "VIS005" - pdict["vis_06"] = "VIS006" - pdict["vis_08"] = "VIS008" - pdict["vis_09"] = "VIS009" - pdict["nir_13"] = "IR_013" - pdict["nir_16"] = "IR_016" - pdict["nir_22"] = "IR_022" - pdict["ir_38"] = "IR_038" - pdict["wv_63"] = "WV_063" - pdict["wv_73"] = "WV_073" - pdict["ir_87"] = "IR_087" - pdict["ir_97"] = "IR_097" - pdict["ir_105"] = "IR_105" - pdict["ir_123"] = "IR_123" - pdict["ir_133"] = "IR_133" - - # For msg, the satpy and icare channel names are the same. - pdict["VIS006"] = "VIS006" - pdict["VIS008"] = "VIS008" - pdict["HRV"] = "HRV" - pdict["IR_016"] = "IR_016" - pdict["IR_039"] = "IR_039" - pdict["WV_062"] = "WV_062" - pdict["WV_073"] = "WV_073" - pdict["IR_087"] = "IR_087" - pdict["IR_097"] = "IR_097" - pdict["IR_108"] = "IR_108" - pdict["IR_120"] = "IR_120" - pdict["IR_134"] = "IR_134" - - # For the goesr satellites : - pdict["C01"] = "VIS_004" - pdict["C02"] = "VIS_006" - pdict["C03"] = "VIS_008" - pdict["C04"] = "VIS_014" - pdict["C05"] = "VIS_016" - pdict["C06"] = "VIS_022" - pdict["C07"] = "IR_039" - pdict["C08"] = "IR_062" - pdict["C09"] = "IR_069" - pdict["C10"] = "IR_073" - pdict["C11"] = "IR_085" - pdict["C12"] = "IR_096" - pdict["C13"] = "IR_103" - pdict["C14"] = "IR_114" - pdict["C15"] = "IR_123" - pdict["C16"] = "IR_133" - - # For himawari. - # BO1 : name in satpy. VIS004 : name in icare/meteofrance netcdf. - pdict["B01"] = "VIS004" - pdict["B02"] = "VIS005" - pdict["B03"] = "VIS006" - pdict["B04"] = "VIS008" - pdict["B05"] = "IR_016" - pdict["B06"] = "IR_022" - pdict["B07"] = "IR_038" - pdict["B08"] = "WV_062" - pdict["B09"] = "WV_069" - pdict["B10"] = "WV_073" - pdict["B11"] = "IR_085" - pdict["B12"] = "IR_096" - pdict["B13"] = "IR_104" - pdict["B14"] = "IR_112" - pdict["B15"] = "IR_123" - pdict["B16"] = "IR_132" - return pdict - # buildChannelCorrespondanceName() - - def _get_dsname(self, ds_id) : - """Return the correct dataset name based on requested band. - ds_id = DataID(name='vis_08', - wavelength=WavelengthRange(...), - resolution=2000, calibration=, - modifiers=()) - """ - pdict = self.buildChannelCorrespondanceName() - - satpyName = ds_id["name"] - if satpyName in pdict : - icareName = pdict[satpyName] - else : - print( - "Soft not adaptated for this channel : ds_id = ", - satpyName) - exit(1) - - return icareName - # _get_dsname() - - def initialisation_dataset(self) : - listeToutesVariables = { - "VIS004", "VIS_004", "VIS005", "VIS_005", "VIS006", - "VIS_006", "VIS006", "VIS_008", "VIS008", "VIS_009", - "VIS009", "IR_013", "IR_016", "IR_022", "IR_038", - "IR_039", "IR_062", "WV_062", "WV_063", "IR_069", - "WV_069", "IR_073", "WV_073", "IR_085", "IR_087", - "IR_096", "IR_097", "IR_103", "IR_104", "IR_105", - "IR_108", "IR_112", "IR_120", "IR_123", "IR_132", - "IR_133", "IR_134", "HRV" - } - - self.mda = {} - self.scale_factor = {} - self.offset = {} - self.alpha = {} - self.beta = {} - self.nuc = {} - self.bandfactor = {} - self.variableComptes = {} - - # Loop over the all possible ds_get_name of every satellites : - for ds_get_name in listeToutesVariables : - if ds_get_name in self.nc : - - variable = self.nc[ds_get_name] - attributs = variable.attrs - - self.scale_factor[ds_get_name] = attributs["scale_factor"] - self.offset[ds_get_name] = attributs["add_offset"] - - if "nuc" in attributs : - # Brightness temperature. - self.alpha[ds_get_name] = attributs["alpha"] - self.beta[ds_get_name] = attributs["beta"] - self.nuc[ds_get_name] = attributs["nuc"] - - nomRetourComptes = "Temp_to_Native_count_" + ds_get_name - - elif "bandfactor" in attributs : - # Albedo. - self.bandfactor[ds_get_name] = attributs["bandfactor"] - - nomRetourComptes = "Albedo_to_Native_count_" + ds_get_name - - else : - print( - "Nuc or bandfactor not founded int the attributs of ", - ds_get_name) - exit(1) - - self.variableComptes[ds_get_name] = self.nc[nomRetourComptes] - # (65536). Correspondence from 0 to 65535 towards - # the original spatial agency counts. - - self.mda[ds_get_name] = self.preparer_metadata(variable) - # initialisation_dataset() - - def comebacktoNativeData(self, ds_get_name) : - """ Come back to the original counts of the hrit. - ds_get_name : meteo france name of a channel : IR_108. """ - - variable = self.nc[ds_get_name] - # Variable is between -9000 to 4000 (temperature) - # or between 0 to 10000 (albedo). - - offset = self.offset[ds_get_name] # 0 or 273.15 - variable += 32768 # 0 to 65535 - - if offset == 0. : - # Albedo. - nom = "Albedo_to_Native_count_" + ds_get_name - else : - nom = "Temp_to_Native_count_" + ds_get_name - """ Temp_to_Native_count_IR_062 """ - - variableComptes = self.nc[nom] # (65536). - # Correspondence from 0 to 65535 towards the original spatial agency counts. - - arrayTableConversion = xr.DataArray.to_numpy(variableComptes) - - tableau = arrayTableConversion[variable[:]] - """ Come back to the original counts of the hrit. - tableau : 0 to 4095 if native datas coded with 12 bits. """ - - variable[:] = tableau - return(variable) - # comebacktoNativeData(self, ds_get_name) - - def comebacktoRadiance(self, ds_get_name) : - # Come back to the radiance. - - variable = self.nc[ds_get_name] - # Variable is between -9000 to 4000 (temperature) - # or between 0 to 10000 (albedo). - - scale_factor = self.scale_factor[ds_get_name] # 0.01 - offset = self.offset[ds_get_name] # 0 or 273.15 - - if offset == 0. : - # Visible channel. - bandfactor = self.bandfactor[ds_get_name] - - # Variable is an albedo from 0 to 10000. - variable = variable * scale_factor / 100. * bandfactor - # => variable is a reflectance between 0 and 1. - # radiance in mWm-2sr-1(cm-1)-1 - else : - # Brightness temperature. - nuc = self.nuc[ds_get_name] - alpha = self.alpha[ds_get_name] - beta = self.beta[ds_get_name] - - variable = variable * scale_factor + offset - # variable is in Kelvin. - variable = variable * alpha + beta - - c1 = 1.1910427e-5 # Planck - c2 = 1.4387752 # Planck - resul1 = c1 * np.power(nuc, 3.) - resul2 = c2 * nuc - val2 = np.exp(resul2 / variable) - 1. - variable = resul1 / val2 - # Radiance in mWm-2sr-1(cm-1)-1 - return(variable) - # comebacktoRadiance(self, ds_get_name) - - def get_dataset(self, ds_id, ds_info) : - """Get the dataset. - ds_id["calibration"] = key["calibration"] = - ["brightness_temperature", "reflectance", "radiance", "counts"] - """ - ds_get_name = self._get_dsname(ds_id) - # "IR_096" - - mda = self.mda[ds_get_name] - mda.update(ds_info) - - calibration = ds_id["calibration"] - - if calibration == "counts" : - # Come back to the original counts of the hrit... - variable = self.comebacktoNativeData(ds_get_name) - - elif calibration == "radiance" : - # Come back to the radiance. - variable = self.comebacktoRadiance(ds_get_name) - - elif calibration == "brightness_temperature" : - variable = self.nc[ds_get_name] - # WV_062 calibration.brightness_temperature, from -9000 to 4000 - scale_factor = self.scale_factor[ds_get_name] - offset = self.offset[ds_get_name] - if offset != 273.15 : - print(ds_get_name, ". offset = ", offset) - message = "Soft not intended for a reflectance " - message += "with a wave length more than 3 microns." - print(message) - exit(1) - - variable = variable * scale_factor + offset - # variable is in Kelvin. - - elif calibration == "reflectance" : - variable = self.nc[ds_get_name] - # VIS006 calibration.reflectance, from 0 to 10000 - scale_factor = self.scale_factor[ds_get_name] - offset = self.offset[ds_get_name] - if offset != 0. : - print(ds_get_name, ". offset = ", offset) - message = "Soft not intended " - message += "for a brightness temperature " - message += "with a wave length less than 3 microns." - print(message) - exit(1) - - variable = variable * scale_factor - # variable is an albedo between 0 and 100. - - else : - print("Calibration mode not expected : ", calibration) - exit(1) - - variable = variable.rename({variable.dims[0] : "y", variable.dims[1] : "x"}) - variable.attrs.update(mda) - return variable - # get_dataset() - - def orbital_param(self) : - orb_param_dict = { - "orbital_parameters": { - "satellite_actual_longitude": self.longitudeReelle, - "satellite_actual_latitude": 0., - "satellite_actual_altitude": self.alt, - "satellite_nominal_longitude": self.longitudedeprojection, - "satellite_nominal_latitude": 0, - "satellite_nominal_altitude": self.alt, - "projection_longitude": self.longitudedeprojection, - "projection_latitude": 0., - "projection_altitude": self.alt, - } - } - - return orb_param_dict - - def channelType(self, pdict, pdictResoAdesc, pdictResoPid, satellite) : - strNbpix = str(self.nbpix) - if strNbpix in pdictResoAdesc : - pdict["a_desc"] = pdictResoAdesc[strNbpix] - pdict["p_id"] = pdictResoPid[strNbpix] - else : - print( - "ERROR : not expected resolution for ", - satellite, " : ", self.nbpix) - exit(1) - return(pdict) - - def get_area_def(self, ds_id) : - """Get the area def.""" - - pdict = {} - pdict["cfac"] = np.int32(self.cfac) - pdict["lfac"] = np.int32(self.lfac) - pdict["coff"] = np.float32(self.coff) - pdict["loff"] = np.float32(self.loff) - - pdict["a"] = 6378169 - pdict["b"] = 6356583.8 - pdict["h"] = self.alt - pdict["a"] - # 36000000 mètres. - pdict["ssp_lon"] = self.longitudedeprojection - pdict["ncols"] = self.nblig - pdict["nlines"] = self.nbpix - pdict["sweep"] = "y" - - # Force scandir to SEVIRI default, not known from file - pdict["scandir"] = "S2N" - pdict["a_name"] = "geosmsg" - - pdictResoAdesc = {} - pdictResoPid = {} - - if self.sensor == "seviri" : - # msg. - pdict["scandir"] = "N2S" - pdict["a_name"] = "geosmsg" - - pdictResoAdesc["3712"] = "MSG/SEVIRI low resolution channel area" - pdictResoPid["3712"] = "msg_lowres" - pdictResoAdesc["11136"] = "MSG/SEVIRI HRV channel area" - pdictResoPid["11136"] = "msg_hires" - - pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "msg") - - elif self.sensor == "fci" : - # mtg. - pdict["scandir"] = "N2S" - pdict["a_name"] = "geosmtg" - - pdictResoAdesc["5568"] = "MTG 2km channel area" - pdictResoPid["5568"] = "mtg_lowres" - pdictResoAdesc["11136"] = "MTG 1km channel area" - pdictResoPid["11136"] = "mtg_midres" - pdictResoAdesc["22272"] = "MTG 500m channel area" - pdictResoPid["22272"] = "mtg_hires" - - pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "mtg") - - elif self.sensor == "ahi" : - # Himawari. - pdict["scandir"] = "N2S" - pdict["a_name"] = "geoshima" - - pdictResoAdesc["5500"] = "HIMA 2km channel area" - pdictResoPid["5500"] = "hima_lowres" - pdictResoAdesc["11000"] = "HIMA 1km channel area" - pdictResoPid["11000"] = "hima_midres" - pdictResoAdesc["22000"] = "HIMA 500m channel area" - pdictResoPid["22000"] = "hima_hires" - - pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "hima") - - elif self.sensor == "abi" : - # Goesr. - pdict["scandir"] = "N2S" - pdict["a_name"] = "geosgoesr" - pdict["sweep"] = "x" - - pdictResoAdesc["5424"] = "GOESR 2km channel area" - pdictResoPid["5424"] = "goesr_lowres" - pdictResoAdesc["10848"] = "GOESR 1km channel area" - pdictResoPid["10848"] = "goesr_midres" - pdictResoAdesc["21696"] = "GOESR 500m channel area" - pdictResoPid["21696"] = "goesr_hires" - - pdict = self.channelType(pdict, pdictResoAdesc, pdictResoPid, "goesr") - - else : - print("ERROR : " + self.sensor + " not expected.") - exit(1) - - aex = get_area_extent(pdict) - area = get_area_definition(pdict, aex) - - # print("area = ", area) - return area - # get_area_def() + # Cf readers/file_handlers.py. + + def __init__(self, filename, filename_info, filetype_info) : + """Init the file handler.""" + + super().__init__(filename, filename_info, filetype_info) + + self.nc = xr.open_dataset( + self.filename, decode_cf=True, mask_and_scale=False, + chunks={"xc": "auto", "yc": "auto"}) + self.metadata = {} + + self.metadata["start_time"] = self.get_endOrStartTime( + "time_coverage_start") + self.metadata["end_time"] = self.get_endOrStartTime( + "time_coverage_end") + + # message = "Reading: " + filename + # message += " start: " + format(self.start_time) + # message += " end: " + format(self.end_time) + # logger.info(message) + + self.netcdfCommonAttributReading() + # __init__() + + def netcdfCommonAttributReading(self) : + self.sensor = self.sensor_name() + # seviri + + self.platform = self.platform_name() + # Meteosat-10 + + self.res() + # Resolution : 3000.4 m + + self.actualLongitude = self.satlon() + self.projectionLongitude = self.projlon() + + self.zone = self.nc.attrs["Area_of_acquisition"] + # globe. + + def sensor_name(self) : + """Get the sensor name seviri, fci, abi, ahi. + """ + variable = self.nc["satellite"] + platform = variable.attrs["id"] # msg1, msg01, MSG1... + platform = platform[:3] # msg, MSG + platform = platform.lower() # msg + + pdict = {} + pdict["msg"] = "seviri" + pdict["mtg"] = "fci" + pdict["goe"] = "abi" + pdict["him"] = "ahi" + + if platform in pdict : + sensor = pdict[platform] + else : + message = "Unsupported satellite platform : " + message += self.platform + raise NotImplementedError(message) + return sensor + + def platform_name(self) : + # Icare and météo france use non-standard platform names. + # Change is needed for pyspectral : + # pyspectral/rsr_seviri_Meteosat-10.h5 in the call + # Calculator(platform_name, sensor, name). + + variable = self.nc["satellite"] + platform = variable.attrs["id"] # msg1, msg01, MSG1... + + pdict = {} + pdict["msg1"] = "Meteosat-08" + pdict["msg01"] = "Meteosat-08" + pdict["MSG1"] = "Meteosat-08" + pdict["msg2"] = "Meteosat-09" + pdict["msg02"] = "Meteosat-09" + pdict["MSG2"] = "Meteosat-09" + pdict["msg3"] = "Meteosat-10" + pdict["msg03"] = "Meteosat-10" + pdict["MSG3"] = "Meteosat-10" + pdict["msg4"] = "Meteosat-11" + pdict["msg04"] = "Meteosat-11" + pdict["MSG4"] = "Meteosat-11" + pdict["mtgi1"] = "Meteosat-12" + pdict["mtg1"] = "Meteosat-12" + pdict["MTG01"] = "Meteosat-12" + pdict["goes16"] = "GOES-16" + pdict["goes17"] = "GOES-17" + pdict["goes18"] = "GOES-18" + pdict["goes19"] = "GOES-19" + pdict["hima08"] = "Himawari-8" + pdict["hima09"] = "Himawari-9" + + if platform in pdict : + platform = pdict[platform] + else : + message = "Unsupported satellite platform : " + platform + raise NotImplementedError(message) + return platform + # platform_name() + + def satlon(self) : + """Get the satellite longitude.""" + variable = self.nc["satellite"] + actualLongitude = variable.attrs["lon"] + return actualLongitude + + def projlon(self): + """Get the projection longitude.""" + variable = self.nc["geos"] + projectionLongitude = variable.attrs["longitude_of_projection_origin"] + return projectionLongitude + + @property + def projection(self): + """Get the projection.""" + return "geos" + + def res(self) : + """Get the resolution. + The resolution can be read in the attribute geotransform + of the following variables : + GeosCoordinateSystem500m, GeosCoordinateSystem_h, + GeosCoordinateSystem1km, GeosCoordinateSystem2km, + GeosCoordinateSystem. + cfac, lfac, coff, loff can be read in the variables ImageNavigationxxx. + """ + + if "GeosCoordinateSystem500m" in self.nc : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem500m"] + Xvariable = self.nc["X500m"] + Yvariable = self.nc["Y500m"] + navigationString = "ImageNavigation500m" + + elif "GeosCoordinateSystem_h" in self.nc : + # Hrv from msg. + variable = self.nc["GeosCoordinateSystem_h"] + Xvariable = self.nc["X_h"] + Yvariable = self.nc["Y_h"] + navigationString = "ImageNavigation_h" + + elif "GeosCoordinateSystem1km" in self.nc : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem1km"] + Xvariable = self.nc["X1km"] + Yvariable = self.nc["Y1km"] + navigationString = "ImageNavigation1km" + + elif "GeosCoordinateSystem2km" in self.nc : + # Mtg, himawari, goesr. + variable = self.nc["GeosCoordinateSystem2km"] + Xvariable = self.nc["X2km"] + Yvariable = self.nc["Y2km"] + navigationString = "ImageNavigation2km" + + elif "GeosCoordinateSystem" in self.nc : + # Msg in 3 kms. + variable = self.nc["GeosCoordinateSystem"] + Xvariable = self.nc["X"] + Yvariable = self.nc["Y"] + navigationString = "ImageNavigation" + + else : + message = "Variables GeosCoordinateSystemXX not founded." + raise NotImplementedError(message) + + geotransform = variable.attrs["GeoTransform"] + # geotransform = -5570254, 3000.40604, 0, 5570254, 0, -3000.40604 + + chunksGeotransform = geotransform.split(", ") + self.resolution = float(chunksGeotransform[1]) + # 3000.40604 + + self.X = Xvariable[:] + self.nbpix = self.X.shape[0] + self.Y = Yvariable[:] + self.nblig = self.Y.shape[0] + + variable = self.nc[navigationString] + self.cfac = float(variable.attrs["CFAC"]) + self.lfac = float(variable.attrs["LFAC"]) + self.coff = float(variable.attrs["COFF"]) + self.loff = float(variable.attrs["LOFF"]) + # res() + + def get_endOrStartTime(self, AttributeName) : + """Get the end or the start time. Global attribute of the netcdf. + AttributName : "time_coverage_start", "time_coverage_end" + """ + attr = self.nc.attrs[AttributeName] + # YYYY-MM-DDTHH:MM:SSZNNN or YYYY-MM-DDTHH:MM:SSZ + # In some versions milliseconds are present, sometimes not. + lengthString = len(attr) + if lengthString == 22 : + # Goesr : 2024-06-28T10:00:21.1Z + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%S.%fZ") + elif lengthString == 20 : + # Mtg. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ") + else : + # Msg, hima. + stacq = dt.datetime.strptime(attr, "%Y-%m-%dT%H:%M:%SZ%f") + return stacq + # get_endOrStartTime() + + @property + def start_time(self) : + return(self.metadata["start_time"]) + + @property + def end_time(self) : + return(self.metadata["end_time"]) + + @property + def alt(self) : + """Get the altitude.""" + variable = self.nc["satellite"] + altitude = variable.attrs["dst"] # 36000000 meters. + altitude += 6378169. # equatorial radius of the earth. + return altitude + + def prepare_metadata(self, variable) : + """Get the metadata for a channel variable. + Add the global attributes.""" + mda = {} + + attributs = variable.attrs + for name in attributs : + mda.update({name: attributs.get(name)}) + + mda.update({ + "start_time": self.start_time, + "end_time": self.end_time, + "platform_name": self.platform, + "sensor": self.sensor, + "zone": self.zone, + "projection_altitude": self.alt, + "cfac": self.cfac, + "lfac": self.lfac, + "coff": self.coff, + "loff": self.loff, + "resolution": self.resolution, + "satellite_actual_longitude": self.actualLongitude, + "projection_longitude": self.projectionLongitude, + "projection_type": self.projection + }) + + mda.update(self.orbital_param()) + return mda + # prepare_metadata(). + + def buildChannelCorrespondanceName(self) : + pdict = {} + + # For mtg. + # vis_04 is the name in satpy. + # VIS004 is the icare/meteofrance netcdf name. + pdict["vis_04"] = "VIS004" + pdict["vis_05"] = "VIS005" + pdict["vis_06"] = "VIS006" + pdict["vis_08"] = "VIS008" + pdict["vis_09"] = "VIS009" + pdict["nir_13"] = "IR_013" + pdict["nir_16"] = "IR_016" + pdict["nir_22"] = "IR_022" + pdict["ir_38"] = "IR_038" + pdict["wv_63"] = "WV_063" + pdict["wv_73"] = "WV_073" + pdict["ir_87"] = "IR_087" + pdict["ir_97"] = "IR_097" + pdict["ir_105"] = "IR_105" + pdict["ir_123"] = "IR_123" + pdict["ir_133"] = "IR_133" + + # For msg, the satpy and icare channel names are the same. + pdict["VIS006"] = "VIS006" + pdict["VIS008"] = "VIS008" + pdict["HRV"] = "HRV" + pdict["IR_016"] = "IR_016" + pdict["IR_039"] = "IR_039" + pdict["WV_062"] = "WV_062" + pdict["WV_073"] = "WV_073" + pdict["IR_087"] = "IR_087" + pdict["IR_097"] = "IR_097" + pdict["IR_108"] = "IR_108" + pdict["IR_120"] = "IR_120" + pdict["IR_134"] = "IR_134" + + # For the goesr satellites : + pdict["C01"] = "VIS_004" + pdict["C02"] = "VIS_006" + pdict["C03"] = "VIS_008" + pdict["C04"] = "VIS_014" + pdict["C05"] = "VIS_016" + pdict["C06"] = "VIS_022" + pdict["C07"] = "IR_039" + pdict["C08"] = "IR_062" + pdict["C09"] = "IR_069" + pdict["C10"] = "IR_073" + pdict["C11"] = "IR_085" + pdict["C12"] = "IR_096" + pdict["C13"] = "IR_103" + pdict["C14"] = "IR_114" + pdict["C15"] = "IR_123" + pdict["C16"] = "IR_133" + + # For himawari. + # BO1 : name in satpy. VIS004 : name in icare/meteofrance netcdf. + pdict["B01"] = "VIS004" + pdict["B02"] = "VIS005" + pdict["B03"] = "VIS006" + pdict["B04"] = "VIS008" + pdict["B05"] = "IR_016" + pdict["B06"] = "IR_022" + pdict["B07"] = "IR_038" + pdict["B08"] = "WV_062" + pdict["B09"] = "WV_069" + pdict["B10"] = "WV_073" + pdict["B11"] = "IR_085" + pdict["B12"] = "IR_096" + pdict["B13"] = "IR_104" + pdict["B14"] = "IR_112" + pdict["B15"] = "IR_123" + pdict["B16"] = "IR_132" + return pdict + # buildChannelCorrespondanceName() + + def _get_dsname(self, ds_id) : + """Return the correct dataset name based on requested band. + ds_id = DataID(name='vis_08', + wavelength=WavelengthRange(...), + resolution=2000, calibration=, + modifiers=()) + """ + pdict = self.buildChannelCorrespondanceName() + + satpyName = ds_id["name"] + if satpyName in pdict : + icareName = pdict[satpyName] + else : + message = "Soft not adaptated for this channel : ds_id = " + message += satpyName + raise NotImplementedError(message) + + return icareName + # _get_dsname() + + def channelAttributs(self, ds_get_name) : + if ds_get_name not in self.nc : + message = "Channel " + ds_get_name + "not founded " + message += "in the netcdf." + raise NotImplementedError(message) + + self.mda = {} + self.scale_factor = {} + self.offset = {} + self.alpha = {} + self.beta = {} + self.nuc = {} + self.bandfactor = {} + self.backtocountVariable = {} + + variable = self.nc[ds_get_name] + attributs = variable.attrs + + self.scale_factor[ds_get_name] = attributs["scale_factor"] + self.offset[ds_get_name] = attributs["add_offset"] + + if "nuc" in attributs : + # Brightness temperature. + self.alpha[ds_get_name] = attributs["alpha"] + self.beta[ds_get_name] = attributs["beta"] + self.nuc[ds_get_name] = attributs["nuc"] + + backtocountName = "Temp_to_Native_count_" + ds_get_name + + elif "bandfactor" in attributs : + # Albedo. + self.bandfactor[ds_get_name] = attributs["bandfactor"] + backtocountName = "Albedo_to_Native_count_" + ds_get_name + + else : + message = "Nuc or bandfactor not founded in the attributs" + message += " of " + ds_get_name + raise NotImplementedError(message) + + self.backtocountVariable[ds_get_name] = self.nc[backtocountName] + # (65536). Correspondence from 0 to 65535 towards + # the original spatial agency counts. + + self.mda[ds_get_name] = self.prepare_metadata(variable) + # channelAttributs(ds_get_name) + + def comebacktoNativeData(self, ds_get_name) : + """ Come back to the original counts of the hrit. + ds_get_name : meteo france name of a channel : IR_108. """ + + variable = self.nc[ds_get_name] + # Variable is between -9000 to 4000 (temperature) + # or between 0 to 10000 (albedo). + + offset = self.offset[ds_get_name] # 0 or 273.15 + variable += 32768 # 0 to 65535 + + if offset == 0. : + # Albedo. + name = "Albedo_to_Native_count_" + ds_get_name + else : + name = "Temp_to_Native_count_" + ds_get_name + """ Temp_to_Native_count_IR_062 """ + + backtocountVariable = self.nc[name] # (65536). + # Correspondence from 0 to 65535 + # towards the original spatial agency counts. + + arrayTableConversion = xr.DataArray.to_numpy(backtocountVariable) + + tableContents = arrayTableConversion[variable[:]] + """ Come back to the original counts of the hrit. + tableau : 0 to 4095 if native datas coded with 12 bits. """ + + variable[:] = tableContents + return(variable) + # comebacktoNativeData(self, ds_get_name) + + def comebacktoRadiance(self, ds_get_name) : + """ Come back to the radiance. + ds_get_name : meteo france name of a channel : IR_108. """ + + variable = self.nc[ds_get_name] + # Variable is between -9000 to 4000 (temperature) + # or between 0 to 10000 (albedo). + + scale_factor = self.scale_factor[ds_get_name] # 0.01 + offset = self.offset[ds_get_name] # 0 or 273.15 + + if offset == 0. : + # Visible channel. + bandfactor = self.bandfactor[ds_get_name] + + # Variable is an albedo from 0 to 10000. + variable = variable * scale_factor / 100. * bandfactor + # => variable is a reflectance between 0 and 1. + # radiance in mWm-2sr-1(cm-1)-1 + else : + # Brightness temperature. + nuc = self.nuc[ds_get_name] + alpha = self.alpha[ds_get_name] + beta = self.beta[ds_get_name] + + variable = variable * scale_factor + offset + # variable becomes Kelvin. + + variable = variable * alpha + beta + resul1 = C1 * np.power(nuc, 3.) + resul2 = C2 * nuc + val2 = np.exp(resul2 / variable) - 1. + variable = resul1 / val2 + # Radiance in mWm-2sr-1(cm-1)-1 + return(variable) + # comebacktoRadiance(self, ds_get_name) + + def get_dataset(self, ds_id, ds_info) : + """Get the dataset. + ds_id["calibration"] = key["calibration"] = + ["brightness_temperature", "reflectance", "radiance", "counts"] + """ + ds_get_name = self._get_dsname(ds_id) + # ds_get_name is the meteo France Icare name of the channel : IR_096. + + self.channelAttributs(ds_get_name) + + mda = self.mda[ds_get_name] + mda.update(ds_info) + + calibration = ds_id["calibration"] + + if calibration == "counts" : + # Come back to the original counts of the hrit... + variable = self.comebacktoNativeData(ds_get_name) + + elif calibration == "radiance" : + # Come back to the radiance. + variable = self.comebacktoRadiance(ds_get_name) + + elif calibration == "brightness_temperature" : + variable = self.nc[ds_get_name] + # WV_062 calibration.brightness_temperature, from -9000 to 4000 + scale_factor = self.scale_factor[ds_get_name] + offset = self.offset[ds_get_name] + if offset != 273.15 : + message = "Soft not intended for a reflectance " + message += "with a wave length more than 3 microns. " + message += ds_get_name + " offset = " + str(offset) + raise NotImplementedError(message) + + variable = variable * scale_factor + offset + # variable becomes Kelvin. + + elif calibration == "reflectance" : + variable = self.nc[ds_get_name] + # VIS006 calibration.reflectance, from 0 to 10000 + scale_factor = self.scale_factor[ds_get_name] + offset = self.offset[ds_get_name] + if offset != 0. : + message = "Soft not intended " + message += "for a brightness temperature " + message += "with a wave length less than 3 microns. " + message += ds_get_name + " offset = " + str(offset) + raise NotImplementedError(message) + + variable = variable * scale_factor + # variable becomes an albedo between 0 and 100. + + else : + message = "Calibration mode not expected : " + calibration + raise NotImplementedError(message) + + variable = variable.rename( + {variable.dims[0] : "y", variable.dims[1] : "x"}) + variable.attrs.update(mda) + return variable + # get_dataset() + + def orbital_param(self) : + orb_param_dict = { + "orbital_parameters": { + "satellite_actual_longitude": self.actualLongitude, + "satellite_actual_latitude": 0., + "satellite_actual_altitude": self.alt, + "satellite_nominal_longitude": self.projectionLongitude, + "satellite_nominal_latitude": 0, + "satellite_nominal_altitude": self.alt, + "projection_longitude": self.projectionLongitude, + "projection_latitude": 0., + "projection_altitude": self.alt, + } + } + return orb_param_dict + + def channelType(self, pdict, pdictResoAdesc, pdictResoPid, satellite) : + strNbpix = str(self.nbpix) + if strNbpix in pdictResoAdesc : + pdict["a_desc"] = pdictResoAdesc[strNbpix] + pdict["p_id"] = pdictResoPid[strNbpix] + else : + message = "Resolution " + str(self.nbpix) + message += " not expected for " + satellite + raise NotImplementedError(message) + + return(pdict) + + def get_area_def(self, ds_id) : + """Get the area def.""" + + pdict = {} + pdict["cfac"] = np.int32(self.cfac) + pdict["lfac"] = np.int32(self.lfac) + pdict["coff"] = np.float32(self.coff) + pdict["loff"] = np.float32(self.loff) + + pdict["a"] = 6378169 + pdict["b"] = 6356583.8 + pdict["h"] = self.alt - pdict["a"] + # 36000000 mètres. + pdict["ssp_lon"] = self.projectionLongitude + pdict["ncols"] = self.nblig + pdict["nlines"] = self.nbpix + pdict["sweep"] = "y" + + # Force scandir to SEVIRI default, not known from file + pdict["scandir"] = "S2N" + pdict["a_name"] = "geosmsg" + + pdictResoAdesc = {} + pdictResoPid = {} + + if self.sensor == "seviri" : + # msg. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geosmsg" + + pdictResoAdesc["3712"] = "MSG/SEVIRI low resolution channel area" + pdictResoPid["3712"] = "msg_lowres" + pdictResoAdesc["11136"] = "MSG/SEVIRI HRV channel area" + pdictResoPid["11136"] = "msg_hires" + + pdict = self.channelType( + pdict, pdictResoAdesc, pdictResoPid, "msg") + + elif self.sensor == "fci" : + # mtg. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geosmtg" + + pdictResoAdesc["5568"] = "MTG 2km channel area" + pdictResoPid["5568"] = "mtg_lowres" + pdictResoAdesc["11136"] = "MTG 1km channel area" + pdictResoPid["11136"] = "mtg_midres" + pdictResoAdesc["22272"] = "MTG 500m channel area" + pdictResoPid["22272"] = "mtg_hires" + + pdict = self.channelType( + pdict, pdictResoAdesc, pdictResoPid, "mtg") + + elif self.sensor == "ahi" : + # Himawari. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geoshima" + + pdictResoAdesc["5500"] = "HIMA 2km channel area" + pdictResoPid["5500"] = "hima_lowres" + pdictResoAdesc["11000"] = "HIMA 1km channel area" + pdictResoPid["11000"] = "hima_midres" + pdictResoAdesc["22000"] = "HIMA 500m channel area" + pdictResoPid["22000"] = "hima_hires" + + pdict = self.channelType( + pdict, pdictResoAdesc, pdictResoPid, "hima") + + elif self.sensor == "abi" : + # Goesr. + pdict["scandir"] = "N2S" + pdict["a_name"] = "geosgoesr" + pdict["sweep"] = "x" + + pdictResoAdesc["5424"] = "GOESR 2km channel area" + pdictResoPid["5424"] = "goesr_lowres" + pdictResoAdesc["10848"] = "GOESR 1km channel area" + pdictResoPid["10848"] = "goesr_midres" + pdictResoAdesc["21696"] = "GOESR 500m channel area" + pdictResoPid["21696"] = "goesr_hires" + + pdict = self.channelType( + pdict, pdictResoAdesc, pdictResoPid, "goesr") + + else : + message = "Sensor " + self.sensor + " not expected." + raise NotImplementedError(message) + + aex = get_area_extent(pdict) + area = get_area_definition(pdict, aex) + + return area + # get_area_def() diff --git a/satpy/tests/reader_tests/test_geos_netcdficare.py b/satpy/tests/reader_tests/test_geos_netcdficare.py index 75e27e6dd7..27b95f900d 100644 --- a/satpy/tests/reader_tests/test_geos_netcdficare.py +++ b/satpy/tests/reader_tests/test_geos_netcdficare.py @@ -17,22 +17,22 @@ # satpy. If not, see . """Tests for the Icare MeteoFrance netcdfs reader, - satpy/readers/geos_netcdficare.py. + satpy/readers/geos_netcdficare.py. SYNTAXE : - pytest + pytest - Before, copy the satpy/tests/reader_tests/test_geos_netcdfcare.py file - into the pytest directory pytest_projec. + Before, copy the satpy/tests/reader_tests/test_geos_netcdfcare.py file + into the pytest directory pytest_projec. EXECUTION TIME : - 20 seconds. + 20 seconds. DATE OF CREATION : - 2024 11th october. + 2024 16th october. LAST VERSIONS : AUTHOR : - Meteo France. + Meteo France. """ import os @@ -45,270 +45,274 @@ from datetime import datetime from netCDF4 import Dataset +import logging +logger = logging.getLogger('netcdficare') + class TestGeosNetcdfIcareReader() : - # Test of the geos_netcdficare reader. - # This reader has been build for the Icare Meteo France netcdfs. - - def test_geosNetcdficare(self, tmp_path) : - """ A dummy netcdf is built. - A scene self.scn for the convection product for the same date - is built. We check that the scene parameters are the same - as thoses in the dummy netcdf. - This procedure is called by pytest. - """ - - self.init(tmp_path) - - startTime = self.scn.start_time - startTimeString = startTime.strftime('%Y-%m-%dT%H:%M:%S') - # 2024-06-28T10:00:40 - assert startTimeString == self.expectedStartTime - - endTime = self.scn.end_time - endTimeString = endTime.strftime('%Y-%m-%dT%H:%M:%S') - # 2024-06-28T10:12:41 - assert endTimeString == self.expectedEndTime - - sensor = self.scn.sensor_names - for isensor in sensor : - capteur = isensor - assert capteur == self.expectedSensor - - platform = "erreur" - altitude = -1. - longitude = 999. - - for data_arr in self.values : - # values come from the scene. - # print("data_arr.attrs = ", data_arr.attrs) - if "platform_name" in data_arr.attrs : - platform = data_arr.attrs["platform_name"] - if "orbital_parameters" in data_arr.attrs : - subAttr = data_arr.attrs["orbital_parameters"] - if "satellite_actual_altitude" in subAttr : - altitude = subAttr["satellite_actual_altitude"] - if "satellite_actual_longitude" in data_arr.attrs : - longitude = data_arr.attrs["satellite_actual_longitude"] - - longitude = float(int(longitude * 10.)) / 10. - assert platform == self.expectedPlatform - assert longitude == self.expectedLongitude - assert altitude == self.expectedAltitude - - xr = self.scn.to_xarray_dataset() - # print("xr = ", xr) - matrice = xr["convection"] - nblin = matrice.shape[1] - nbpix = matrice.shape[2] - assert nblin == self.expectedNblin - assert nbpix == self.expectedNbpix - - cfac = xr.attrs["cfac"] - assert cfac == self.expectedCfac - lfac = xr.attrs["lfac"] - assert lfac == self.expectedLfac - coff = xr.attrs["coff"] - assert coff == self.expectedCoff - loff = xr.attrs["loff"] - assert loff == self.expectedLoff - - satpyId = xr.attrs["_satpy_id"] - # DataID(name='convection', resolution=3000.403165817) - # Cf satpy/dataset/dataid.py. - - resolution = satpyId.get("resolution") - resolution = float(int(resolution * 10.)) / 10. - assert resolution == self.expectedResolution - - # A picture of convection composite will be displayed. - self.scn.show("convection") - print("The picture should be pink.") - # test_geosNetcdficare(self, tmp_path) - - def init(self, tmp_path) : - """ - A fake netcdf is built. - A scene is built with the reader to be tested, applied to this netcdf. - Called by test_geosNetcdficare(). - """ - self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" - self.filepath = tmp_path - - self.buildNetcdf(self.netcdfName) - - # We will check that the parameters written in the dummy netcdf can be read. - self.expectedStartTime = "2024-06-28T10:00:09" - self.expectedEndTime = "2024-06-28T10:12:41" - actualAltitude = 35786691 + 6378169 # 42164860.0 - actualLongitude = 0.1 - self.expectedPlatform = "Meteosat-10" - self.expectedSensor = "seviri" - self.expectedAltitude = actualAltitude - self.expectedLongitude = actualLongitude - self.expectedCfac = 1.3642337E7 - self.expectedLfac = 1.3642337E7 - self.expectedCoff = 1857.0 - self.expectedLoff = 1857.0 - self.expectedResolution = 3000.4 - self.expectedNblin = 3712 - self.expectedNbpix = 3712 - - # To build a scene at the date 20240628_100000, - # a netcdf corresponding to msg_netcdficare - # is looked for in the filepath directory. - yaml_file = 'msg_netcdficare' - myfiles = find_files_and_readers( - base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), - end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) - print("Found myfiles = ", myfiles) - # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} - - self.scn = Scene(filenames=myfiles, reader=yaml_file) - - print(self.scn.available_dataset_names()) - # ['IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', - # 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073'] - - print(self.scn.available_composite_names()) - """ Static compositor {'_satpy_id': DataID(name='_night_background'), - 'standard_name': 'night_background', - 'prerequisites': [], - 'optional_prerequisites': []} - Static compositor {'_satpy_id': DataID(name='_night_background_hires'), - 'standard_name': 'night_background', 'prerequisites': [], - 'optional_prerequisites': []} - ['airmass', 'ash', 'cloud_phase_distinction', 'cloud_phase_distinction_raw', - 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection'... - """ - - self.scn.load(['convection']) - self.values = self.scn.values() - # init() - - def buildNetcdf(self, ncName) : - """ - ncName : tmp_path / Mmultic3kmNC4_msg03_202406281000.nc - A dummy icare Meteo France netcdf is built here. - Called by init(). - """ - if os.path.exists(ncName) : - os.remove(ncName) - ncfileOut = Dataset( - ncName, mode="w", clobber=True, - diskless=False, persist=False, format='NETCDF4') - - ncfileOut.createDimension(u'ny', 3712) - ncfileOut.createDimension(u'nx', 3712) - ncfileOut.createDimension(u'numerical_count', 65536) - ncfileOut.setncattr("time_coverage_start", "2024-06-28T10:00:09Z383") - ncfileOut.setncattr("time_coverage_end", "2024-06-28T10:12:41Z365") - ncfileOut.setncattr("Area_of_acquisition", "globe") - - fill_value = -32768 - var = ncfileOut.createVariable( - "satellite", "c", zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - - var.setncattr("id", "msg03") - var.setncattr("dst", 35786691.) - var.setncattr("lon", float(0.1)) - - var = ncfileOut.createVariable( - "geos", "c", zlib=True, complevel=4, shuffle=True, - fletcher32=False, contiguous=False, chunksizes=None, - endian='native', least_significant_digit=None) - var.setncattr("longitude_of_projection_origin", 0.) - - var = ncfileOut.createVariable( - "GeosCoordinateSystem", "c", zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr( - "GeoTransform", - "-5570254, 3000.40604, 0, 5570254, 0, -3000.40604") - - var = ncfileOut.createVariable( - "ImageNavigation", "c", zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr("CFAC", 1.3642337E7) - var.setncattr("LFAC", 1.3642337E7) - var.setncattr("COFF", 1857.0) - var.setncattr("LOFF", 1857.0) - - var = ncfileOut.createVariable( - "X", 'float32', u'nx', zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - x0 = -5570254. - dx = 3000.40604 - var[:] = np.array(([(x0 + dx * i) for i in range(3712)])) - - var = ncfileOut.createVariable( - "Y", 'float32', u'ny', zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - y0 = 5570254. - dy = -3000.40604 - var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) - - self.visibleChannelsCreation(ncfileOut, fill_value) - self.infrarougeChannelsCreation(ncfileOut, fill_value) - - ncfileOut.close - # buildNetcdf() - - def visibleChannelsCreation(self, ncfileOut, fill_value) : - for channel in {"VIS006", "VIS008", "IR_016"} : - var = ncfileOut.createVariable( - channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', - least_significant_digit=None, fill_value=fill_value) - var[:] = np.array(([[i * 2 for i in range(3712)] for j in range(3712)])) - # Hundredths of albedo between 0 and 10000. - var.setncattr("scale_factor", 0.01) - var.setncattr("add_offset", 0.) - var.setncattr("bandfactor", 20.76) - var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") - - var = ncfileOut.createVariable( - "Albedo_to_Native_count_" + channel, 'short', - 'numerical_count', zlib=True, complevel=4, shuffle=True, - fletcher32=False, contiguous=False, chunksizes=None, - endian='native', least_significant_digit=None, - fill_value=-9999) - var[:] = np.array(([-9999 for i in range(65536)])) - # In order to come back to the native datas on 10, 12 or 16 bits. - - def infrarougeChannelsCreation(self, ncfileOut, fill_value) : - for channel in { - "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", - "IR_108", "IR_120", "IR_134"} : - var = ncfileOut.createVariable( - channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', - least_significant_digit=None, fill_value=fill_value) - var[:] = np.array( - ([[-9000 + j * 4 for i in range(3712)] for j in range(3712)])) - # Hundredths of celcius degrees. - var.setncattr("scale_factor", 0.01) - var.setncattr("add_offset", 273.15) - var.setncattr("nuc", 1600.548) - var.setncattr("alpha", 0.9963) - var.setncattr("beta", 2.185) - var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") - - var = ncfileOut.createVariable( - "Temp_to_Native_count_" + channel, 'short', - 'numerical_count', zlib=True, complevel=4, shuffle=True, - fletcher32=False, contiguous=False, chunksizes=None, - endian='native', least_significant_digit=None, - fill_value=-9999) - var[:] = np.array(([-9999 for i in range(65536)])) - # In order to come back to the native datas on 10, 12 or 16 bits. - - # class TestGeosNetcdfIcareReader. + # Test of the geos_netcdficare reader. + # This reader has been build for the Icare Meteo France netcdfs. + + def test_geosNetcdficare(self, tmp_path) : + """ A dummy netcdf is built. + A scene self.scn for the convection product for the same date + is built. We check that the scene parameters are the same + as thoses in the dummy netcdf. + This procedure is called by pytest. + """ + + self.init(tmp_path) + + startTime = self.scn.start_time + startTimeString = startTime.strftime('%Y-%m-%dT%H:%M:%S') + # 2024-06-28T10:00:40 + assert startTimeString == self.expectedStartTime + + endTime = self.scn.end_time + endTimeString = endTime.strftime('%Y-%m-%dT%H:%M:%S') + # 2024-06-28T10:12:41 + assert endTimeString == self.expectedEndTime + + sensor = self.scn.sensor_names + for isensor in sensor : + capteur = isensor + assert capteur == self.expectedSensor + + platform = "error" + altitude = -1. + longitude = 999. + + for data_arr in self.values : + # values come from the scene. + # print("data_arr.attrs = ", data_arr.attrs) + if "platform_name" in data_arr.attrs : + platform = data_arr.attrs["platform_name"] + if "orbital_parameters" in data_arr.attrs : + subAttr = data_arr.attrs["orbital_parameters"] + if "satellite_actual_altitude" in subAttr : + altitude = subAttr["satellite_actual_altitude"] + if "satellite_actual_longitude" in data_arr.attrs : + longitude = data_arr.attrs["satellite_actual_longitude"] + + longitude = float(int(longitude * 10.)) / 10. + assert platform == self.expectedPlatform + assert longitude == self.expectedLongitude + assert altitude == self.expectedAltitude + + xr = self.scn.to_xarray_dataset() + matrice = xr["convection"] + nblin = matrice.shape[1] + nbpix = matrice.shape[2] + assert nblin == self.expectedNblin + assert nbpix == self.expectedNbpix + + cfac = xr.attrs["cfac"] + assert cfac == self.expectedCfac + lfac = xr.attrs["lfac"] + assert lfac == self.expectedLfac + coff = xr.attrs["coff"] + assert coff == self.expectedCoff + loff = xr.attrs["loff"] + assert loff == self.expectedLoff + + satpyId = xr.attrs["_satpy_id"] + # DataID(name='convection', resolution=3000.403165817) + # Cf satpy/dataset/dataid.py. + + resolution = satpyId.get("resolution") + resolution = float(int(resolution * 10.)) / 10. + assert resolution == self.expectedResolution + + # A picture of convection composite will be displayed. + self.scn.show("convection") + logger.info("The picture should be pink.") + # test_geosNetcdficare(self, tmp_path) + + def init(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_geosNetcdficare(). + """ + self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" + self.filepath = tmp_path + + self.buildNetcdf(self.netcdfName) + + # We will check that the parameters written in the dummy netcdf + # can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:12:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "Meteosat-10" + self.expectedSensor = "seviri" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 1.3642337E7 + self.expectedLfac = 1.3642337E7 + self.expectedCoff = 1857.0 + self.expectedLoff = 1857.0 + self.expectedResolution = 3000.4 + self.expectedNblin = 3712 + self.expectedNbpix = 3712 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for in the filepath directory. + yaml_file = 'msg_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + logger.info("Found myfiles = ", myfiles) + # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + + # logger.info(self.scn.available_dataset_names()) + # ['IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', + # 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073'] + + # logger.info(self.scn.available_composite_names()) + """ Static compositor {'_satpy_id': DataID(name='_night_background'), + 'standard_name': 'night_background', + 'prerequisites': [], + 'optional_prerequisites': []} + Static compositor {'_satpy_id': DataID(name='_night_background_hires'), + 'standard_name': 'night_background', 'prerequisites': [], + 'optional_prerequisites': []} + ['airmass', 'ash', 'cloud_phase_distinction', ... + 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection'... + """ + + self.scn.load(['convection']) + self.values = self.scn.values() + # init() + + def buildNetcdf(self, ncName) : + """ + ncName : tmp_path / Mmultic3kmNC4_msg03_202406281000.nc + A dummy icare Meteo France netcdf is built here. + Called by init(). + """ + if os.path.exists(ncName) : + os.remove(ncName) + ncfileOut = Dataset( + ncName, mode="w", clobber=True, + diskless=False, persist=False, format='NETCDF4') + + ncfileOut.createDimension(u'ny', 3712) + ncfileOut.createDimension(u'nx', 3712) + ncfileOut.createDimension(u'numerical_count', 65536) + ncfileOut.setncattr("time_coverage_start", "2024-06-28T10:00:09Z383") + ncfileOut.setncattr("time_coverage_end", "2024-06-28T10:12:41Z365") + ncfileOut.setncattr("Area_of_acquisition", "globe") + + fill_value = -32768 + var = ncfileOut.createVariable( + "satellite", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + + var.setncattr("id", "msg03") + var.setncattr("dst", 35786691.) + var.setncattr("lon", float(0.1)) + + var = ncfileOut.createVariable( + "geos", "c", zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None) + var.setncattr("longitude_of_projection_origin", 0.) + + var = ncfileOut.createVariable( + "GeosCoordinateSystem", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + var.setncattr( + "GeoTransform", + "-5570254, 3000.40604, 0, 5570254, 0, -3000.40604") + + var = ncfileOut.createVariable( + "ImageNavigation", "c", zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + var.setncattr("CFAC", 1.3642337E7) + var.setncattr("LFAC", 1.3642337E7) + var.setncattr("COFF", 1857.0) + var.setncattr("LOFF", 1857.0) + + var = ncfileOut.createVariable( + "X", 'float32', u'nx', zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + x0 = -5570254. + dx = 3000.40604 + var[:] = np.array(([(x0 + dx * i) for i in range(3712)])) + + var = ncfileOut.createVariable( + "Y", 'float32', u'ny', zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', least_significant_digit=None) + y0 = 5570254. + dy = -3000.40604 + var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) + + self.visibleChannelsCreation(ncfileOut, fill_value) + self.infrarougeChannelsCreation(ncfileOut, fill_value) + + ncfileOut.close + # buildNetcdf() + + def visibleChannelsCreation(self, ncfileOut, fill_value) : + for channel in {"VIS006", "VIS008", "IR_016"} : + var = ncfileOut.createVariable( + channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', + least_significant_digit=None, fill_value=fill_value) + var[:] = np.array( + ([[i * 2 for i in range(3712)] for j in range(3712)])) + # Hundredths of albedo between 0 and 10000. + var.setncattr("scale_factor", 0.01) + var.setncattr("add_offset", 0.) + var.setncattr("bandfactor", 20.76) + var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + + var = ncfileOut.createVariable( + "Albedo_to_Native_count_" + channel, 'short', + 'numerical_count', zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None, + fill_value=-9999) + var[:] = np.array(([-9999 for i in range(65536)])) + # In order to come back to the native datas on 10, 12 or 16 bits. + + def infrarougeChannelsCreation(self, ncfileOut, fill_value) : + for channel in { + "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", + "IR_108", "IR_120", "IR_134"} : + var = ncfileOut.createVariable( + channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, + shuffle=True, fletcher32=False, contiguous=False, + chunksizes=None, endian='native', + least_significant_digit=None, fill_value=fill_value) + var[:] = np.array( + ([[-9000 + j * 4 for i in range(3712)] for j in range(3712)])) + # Hundredths of celcius degrees. + var.setncattr("scale_factor", 0.01) + var.setncattr("add_offset", 273.15) + var.setncattr("nuc", 1600.548) + var.setncattr("alpha", 0.9963) + var.setncattr("beta", 2.185) + var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + + var = ncfileOut.createVariable( + "Temp_to_Native_count_" + channel, 'short', + 'numerical_count', zlib=True, complevel=4, shuffle=True, + fletcher32=False, contiguous=False, chunksizes=None, + endian='native', least_significant_digit=None, + fill_value=-9999) + var[:] = np.array(([-9999 for i in range(65536)])) + # In order to come back to the native datas on 10, 12 or 16 bits. + + # class TestGeosNetcdfIcareReader. diff --git a/test_geos_netcdficare.py b/test_geos_netcdficare.py deleted file mode 100644 index cf1f5401ad..0000000000 --- a/test_geos_netcdficare.py +++ /dev/null @@ -1,308 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright (c) 2019 Satpy developers -# -# This file is part of satpy. -# -# satpy is free software: you can redistribute it and/or modify it under the -# terms of the GNU General Public License as published by the Free Software -# Foundation, either version 3 of the License, or (at your option) any later -# version. -# -# satpy is distributed in the hope that it will be useful, but WITHOUT ANY -# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR -# A PARTICULAR PURPOSE. See the GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License along with -# satpy. If not, see . - -"""Tests for the Icare MeteoFrance netcdfs reader, - satpy/readers/geos_netcdficare.py. - -SYNTAXE : - pytest - - Before, copy the satpy/tests/reader_tests/test_geos_netcdfcare.py file - into the pytest directory pytest_projec. - -EXECUTION TIME : - 20 seconds. -DATE OF CREATION : - 2024 11th october. -LAST VERSIONS : - -AUTHOR : - Meteo France. -""" - -import os - -import numpy as np - -from satpy.scene import Scene -from satpy import find_files_and_readers - -from datetime import datetime -from netCDF4 import Dataset - - -class TestGeosNetcdfIcareReader() : - # Test of the geos_netcdficare reader. - # This reader has been build for the Icare Meteo France netcdfs. - - def test_geos_netcdficare(self, tmp_path) : - """ A dummy netcdf is built. - A scene self.scn for the convection product for the same date - is built. We check that the scene parameters are the same - as thoses in the dummy netcdf. - This procedure is called by pytest. - """ - - self.init(tmp_path) - - startTime = self.scn.start_time - startTimeString = startTime.strftime('%Y-%m-%dT%H:%M:%S') - # 2024-06-28T10:00:40 - assert startTimeString == self.expectedStartTime - - endTime = self.scn.end_time - endTimeString = endTime.strftime('%Y-%m-%dT%H:%M:%S') - # 2024-06-28T10:12:41 - assert endTimeString == self.expectedEndTime - - sensor = self.scn.sensor_names - for isensor in sensor : - capteur = isensor - assert capteur == self.expectedSensor - - platform = "erreur" - altitude = -1. - longitude = 999. - - for data_arr in self.values : - # values come from the scene. - # print("data_arr.attrs = ", data_arr.attrs) - if "platform_name" in data_arr.attrs : - platform = data_arr.attrs["platform_name"] - if "orbital_parameters" in data_arr.attrs : - subAttr = data_arr.attrs["orbital_parameters"] - if "satellite_actual_altitude" in subAttr : - altitude = subAttr["satellite_actual_altitude"] - if "satellite_actual_longitude" in data_arr.attrs : - longitude = data_arr.attrs["satellite_actual_longitude"] - - longitude = float(int(longitude * 10.)) / 10. - assert platform == self.expectedPlatform - assert longitude == self.expectedLongitude - assert altitude == self.expectedAltitude - - xr = self.scn.to_xarray_dataset() - # print("xr = ", xr) - matrice = xr["convection"] - nblin = matrice.shape[1] - nbpix = matrice.shape[2] - assert nblin == self.expectedNblin - assert nbpix == self.expectedNbpix - - cfac = xr.attrs["cfac"] - assert cfac == self.expectedCfac - lfac = xr.attrs["lfac"] - assert lfac == self.expectedLfac - coff = xr.attrs["coff"] - assert coff == self.expectedCoff - loff = xr.attrs["loff"] - assert loff == self.expectedLoff - - satpyId = xr.attrs["_satpy_id"] - # DataID(name='convection', resolution=3000.403165817) - # Cf satpy/dataset/dataid.py. - - resolution = satpyId.get("resolution") - resolution = float(int(resolution * 10.)) / 10. - assert resolution == self.expectedResolution - - # A picture of convection composite will be displayed. - self.scn.show("convection") - print("The picture should be pink.") - # test_geos_netcdficare(self, tmp_path) - - def init(self, tmp_path) : - """ - A fake netcdf is built. - A scene is built with the reader to be tested, applied to this netcdf. - Called by test_geos_netcdficare(). - """ - self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" - self.filepath = tmp_path - - self.buildNetcdf(self.netcdfName) - - # We will check that the parameters written in the dummy netcdf can be read. - self.expectedStartTime = "2024-06-28T10:00:09" - self.expectedEndTime = "2024-06-28T10:12:41" - actualAltitude = 35786691 + 6378169 # 42164860.0 - actualLongitude = 0.1 - self.expectedPlatform = "Meteosat-10" - self.expectedSensor = "seviri" - self.expectedAltitude = actualAltitude - self.expectedLongitude = actualLongitude - self.expectedCfac = 1.3642337E7 - self.expectedLfac = 1.3642337E7 - self.expectedCoff = 1857.0 - self.expectedLoff = 1857.0 - self.expectedResolution = 3000.4 - self.expectedNblin = 3712 - self.expectedNbpix = 3712 - - # To build a scene at the date 20240628_100000, - # a netcdf corresponding to msg_netcdficare - # is looked for in the filepath directory. - yaml_file = 'msg_netcdficare' - myfiles = find_files_and_readers( - base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), - end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) - print("Found myfiles = ", myfiles) - # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} - - self.scn = Scene(filenames=myfiles, reader=yaml_file) - - print(self.scn.available_dataset_names()) - # ['IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', - # 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073'] - - print(self.scn.available_composite_names()) - """ Static compositor {'_satpy_id': DataID(name='_night_background'), - 'standard_name': 'night_background', - 'prerequisites': [], - 'optional_prerequisites': []} - Static compositor {'_satpy_id': DataID(name='_night_background_hires'), - 'standard_name': 'night_background', 'prerequisites': [], - 'optional_prerequisites': []} - ['airmass', 'ash', 'cloud_phase_distinction', 'cloud_phase_distinction_raw', - 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection'... - """ - - self.scn.load(['convection']) - self.values = self.scn.values() - # init() - - def buildNetcdf(self, ncName) : - """ - ncName : tmp_path / Mmultic3kmNC4_msg03_202406281000.nc - A dummy icare Meteo France netcdf is built here. - Called by init(). - """ - if os.path.exists(ncName) : - os.remove(ncName) - ncfileOut = Dataset( - ncName, mode="w", clobber=True, - diskless=False, persist=False, format='NETCDF4') - - ncfileOut.createDimension(u'ny', 3712) - ncfileOut.createDimension(u'nx', 3712) - ncfileOut.createDimension(u'numerical_count', 65536) - ncfileOut.setncattr("time_coverage_start", "2024-06-28T10:00:09Z383") - ncfileOut.setncattr("time_coverage_end", "2024-06-28T10:12:41Z365") - ncfileOut.setncattr("Area_of_acquisition", "globe") - - fill_value = -32768 - var = ncfileOut.createVariable( - "satellite", "c", zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - - var.setncattr("id", "msg03") - var.setncattr("dst", 35786691.) - var.setncattr("lon", float(0.1)) - - var = ncfileOut.createVariable( - "geos", "c", zlib=True, complevel=4, shuffle=True, - fletcher32=False, contiguous=False, chunksizes=None, - endian='native', least_significant_digit=None) - var.setncattr("longitude_of_projection_origin", 0.) - - var = ncfileOut.createVariable( - "GeosCoordinateSystem", "c", zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr( - "GeoTransform", - "-5570254, 3000.40604, 0, 5570254, 0, -3000.40604") - - var = ncfileOut.createVariable( - "ImageNavigation", "c", zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr("CFAC", 1.3642337E7) - var.setncattr("LFAC", 1.3642337E7) - var.setncattr("COFF", 1857.0) - var.setncattr("LOFF", 1857.0) - - var = ncfileOut.createVariable( - "X", 'float32', u'nx', zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - x0 = -5570254. - dx = 3000.40604 - var[:] = np.array(([(x0 + dx * i) for i in range(3712)])) - - var = ncfileOut.createVariable( - "Y", 'float32', u'ny', zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', least_significant_digit=None) - y0 = 5570254. - dy = -3000.40604 - var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) - - for channel in {"VIS006", "VIS008", "IR_016"} : - var = ncfileOut.createVariable( - channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', - least_significant_digit=None, fill_value=fill_value) - var[:] = np.array(([[i * 2 for i in range(3712)] for j in range(3712)])) - # Hundredths of albedo between 0 and 10000. - var.setncattr("scale_factor", 0.01) - var.setncattr("add_offset", 0.) - var.setncattr("bandfactor", 20.76) - var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") - - var = ncfileOut.createVariable( - "Albedo_to_Native_count_" + channel, 'short', - 'numerical_count', zlib=True, complevel=4, shuffle=True, - fletcher32=False, contiguous=False, chunksizes=None, - endian='native', least_significant_digit=None, - fill_value=-9999) - var[:] = np.array(([-9999 for i in range(65536)])) - # In order to come back to the native datas on 10, 12 or 16 bits. - - for channel in { - "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", - "IR_108", "IR_120", "IR_134"} : - var = ncfileOut.createVariable( - channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, - shuffle=True, fletcher32=False, contiguous=False, - chunksizes=None, endian='native', - least_significant_digit=None, fill_value=fill_value) - var[:] = np.array( - ([[-9000 + j * 4 for i in range(3712)] for j in range(3712)])) - # Hundredths of celcius degrees. - var.setncattr("scale_factor", 0.01) - var.setncattr("add_offset", 273.15) - var.setncattr("nuc", 1600.548) - var.setncattr("alpha", 0.9963) - var.setncattr("beta", 2.185) - var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") - - var = ncfileOut.createVariable( - "Temp_to_Native_count_" + channel, 'short', - 'numerical_count', zlib=True, complevel=4, shuffle=True, - fletcher32=False, contiguous=False, chunksizes=None, - endian='native', least_significant_digit=None, - fill_value=-9999) - var[:] = np.array(([-9999 for i in range(65536)])) - # In order to come back to the native datas on 10, 12 or 16 bits. - ncfileOut.close - # buildNetcdf() - - # class TestGeosNetcdfIcareReader. From 8bebfc138afe6a29bef97963fa82787062beaabc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?laurent=20p=C3=A9rier?= Date: Thu, 28 Nov 2024 13:42:28 +0100 Subject: [PATCH 8/9] Increasing the coverage of tested code. --- .../reader_tests/test_geos_netcdficare.py | 414 +++++++++++++++--- 1 file changed, 342 insertions(+), 72 deletions(-) diff --git a/satpy/tests/reader_tests/test_geos_netcdficare.py b/satpy/tests/reader_tests/test_geos_netcdficare.py index 27b95f900d..992971f44c 100644 --- a/satpy/tests/reader_tests/test_geos_netcdficare.py +++ b/satpy/tests/reader_tests/test_geos_netcdficare.py @@ -26,13 +26,13 @@ into the pytest directory pytest_projec. EXECUTION TIME : - 20 seconds. + 4 minutes. DATE OF CREATION : - 2024 16th october. + 2024 11th october. LAST VERSIONS : AUTHOR : - Meteo France. + Météo France. """ import os @@ -53,16 +53,49 @@ class TestGeosNetcdfIcareReader() : # Test of the geos_netcdficare reader. # This reader has been build for the Icare Meteo France netcdfs. - def test_geosNetcdficare(self, tmp_path) : + def test_mtg_netcdficare(self, tmp_path) : """ A dummy netcdf is built. - A scene self.scn for the convection product for the same date + A scene self.scn for the nir_16 product for the same date is built. We check that the scene parameters are the same as thoses in the dummy netcdf. This procedure is called by pytest. """ + self.initMtg(tmp_path) + self.scn.load(['nir_16']) + self.values = self.scn.values() + self.checkingSceneParameter("nir_16") + # test_mtg_netcdficare() + + def test_msg_netcdficare(self, tmp_path) : + self.initMsgHrv(tmp_path) + self.scn.load(['HRV']) + self.values = self.scn.values() + self.checkingSceneParameter("HRV") + + self.initMsg(tmp_path) + self.scn.load(['convection']) + self.values = self.scn.values() + self.checkingSceneParameter("convection") + # test_msg_netcdficare() - self.init(tmp_path) + def test_hima_netcdficare(self, tmp_path) : + self.initHima(tmp_path) + self.scn.load(['B10']) + self.values = self.scn.values() + self.checkingSceneParameter("B10") + # test_hima_netcdficare() + def test_goesr_netcdficare(self, tmp_path) : + self.initGoesr(tmp_path) + self.scn.load(['C02']) + self.values = self.scn.values() + self.checkingSceneParameter("C02") + # test_goesr_netcdficare() + + # ----------------------------------------------------- # + # typeImage : convection, airmass... + # ----------------------------------------------------- # + def checkingSceneParameter(self, typeImage) : startTime = self.scn.start_time startTimeString = startTime.strftime('%Y-%m-%dT%H:%M:%S') # 2024-06-28T10:00:40 @@ -84,7 +117,6 @@ def test_geosNetcdficare(self, tmp_path) : for data_arr in self.values : # values come from the scene. - # print("data_arr.attrs = ", data_arr.attrs) if "platform_name" in data_arr.attrs : platform = data_arr.attrs["platform_name"] if "orbital_parameters" in data_arr.attrs : @@ -100,9 +132,20 @@ def test_geosNetcdficare(self, tmp_path) : assert altitude == self.expectedAltitude xr = self.scn.to_xarray_dataset() - matrice = xr["convection"] - nblin = matrice.shape[1] - nbpix = matrice.shape[2] + matrice = xr[typeImage] + nbdim = len(matrice.shape) + if nbdim == 3 : + # RGB. + nblin = matrice.shape[1] + nbpix = matrice.shape[2] + elif nbdim == 2 : + # PGM. + nblin = matrice.shape[0] + nbpix = matrice.shape[1] + else : + print("Dimension of shape not expected : ", nbdim) + exit(1) + assert nblin == self.expectedNblin assert nbpix == self.expectedNbpix @@ -122,22 +165,29 @@ def test_geosNetcdficare(self, tmp_path) : resolution = satpyId.get("resolution") resolution = float(int(resolution * 10.)) / 10. assert resolution == self.expectedResolution + # checkingSceneParameter() - # A picture of convection composite will be displayed. - self.scn.show("convection") - logger.info("The picture should be pink.") - # test_geosNetcdficare(self, tmp_path) - - def init(self, tmp_path) : + def initMsg(self, tmp_path) : """ A fake netcdf is built. A scene is built with the reader to be tested, applied to this netcdf. - Called by test_geosNetcdficare(). + Called by test_msg_netcdficare(). """ self.netcdfName = tmp_path / "Mmultic3kmNC4_msg03_202406281000.nc" self.filepath = tmp_path - self.buildNetcdf(self.netcdfName) + listIR = { + "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", + "IR_108", "IR_120", "IR_134" + } + + self.buildNetcdf( + self.netcdfName, 3712, "msg03", x0=-5570254, dx=3000.40604, + y0=5570254, cfac=1.3642337E7, coff=1857.0, + listVisible={"VIS006", "VIS008", "IR_016"}, listIR=listIR, + coordinateSystemName="GeosCoordinateSystem", + nomImageNavigation="ImageNavigation", + nomX="X", nomY="Y", time_coverage_end="2024-06-28T10:12:41Z365") # We will check that the parameters written in the dummy netcdf # can be read. @@ -168,32 +218,248 @@ def init(self, tmp_path) : # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} self.scn = Scene(filenames=myfiles, reader=yaml_file) + # initMsg() - # logger.info(self.scn.available_dataset_names()) - # ['IR_016', 'IR_039', 'IR_087', 'IR_097', 'IR_108', 'IR_120', - # 'IR_134', 'VIS006', 'VIS008', 'WV_062', 'WV_073'] - - # logger.info(self.scn.available_composite_names()) - """ Static compositor {'_satpy_id': DataID(name='_night_background'), - 'standard_name': 'night_background', - 'prerequisites': [], - 'optional_prerequisites': []} - Static compositor {'_satpy_id': DataID(name='_night_background_hires'), - 'standard_name': 'night_background', 'prerequisites': [], - 'optional_prerequisites': []} - ['airmass', 'ash', 'cloud_phase_distinction', ... - 'cloudtop', 'cloudtop_daytime', 'colorized_ir_clouds', 'convection'... + def initMsgHrv(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_msg_netcdficare(). """ + self.netcdfName = tmp_path / "Mmultic1kmNC4_msg03_202406281000.nc" + self.filepath = tmp_path - self.scn.load(['convection']) - self.values = self.scn.values() - # init() + self.buildNetcdf( + self.netcdfName, 11136, "msg03", x0=-5571254., dx=1000.135, + y0=5570254., cfac=40927014, coff=5566.0, listVisible={"HRV"}, + listIR={}, coordinateSystemName="GeosCoordinateSystem_h", + nomImageNavigation="ImageNavigation_h", nomX="X_h", nomY="Y_h", + time_coverage_end="2024-06-28T10:12:41Z365") - def buildNetcdf(self, ncName) : + # We will check that the parameters written in the dummy netcdf + # can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:12:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "Meteosat-10" + self.expectedSensor = "seviri" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 40927014 + self.expectedLfac = 40927014 + self.expectedCoff = 5566.0 + self.expectedLoff = 5566.0 + self.expectedResolution = 1000.1 + self.expectedNblin = 11136 + self.expectedNbpix = 11136 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for, in the filepath directory. + yaml_file = 'msg_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + # logger.info("Found myfiles = ", myfiles) + # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + # initMsgHrv() + + def initHima(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_hima_netcdficare(). + """ + self.netcdfName = tmp_path / "Jmultic2kmNC4_hima09_202406281000.nc" + self.filepath = tmp_path + + listVisible = { + "VIS004", "VIS005", "VIS006", "VIS008", "IR_016", "IR_022" + } + listIR = { + "IR_038", "WV_062", "WV_069", "WV_073", "IR_085", "IR_096", + "IR_104", "IR_112", "IR_123", "IR_132" + } + self.buildNetcdf( + self.netcdfName, 5500, "hima09", x0=-5500000, dx=2000.0000047, + y0=5500000, cfac=20466275., coff=2750.5, + listVisible=listVisible, listIR=listIR, + coordinateSystemName="GeosCoordinateSystem2km", + nomImageNavigation="ImageNavigation2km", + nomX="X2km", nomY="Y2km", + time_coverage_end="2024-06-28T10:08:41Z365") + + # We will check that the parameters written in the dummy netcdf + # can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:08:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "Himawari-9" + self.expectedSensor = "ahi" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 20466275. + self.expectedLfac = 20466275. + self.expectedCoff = 2750.5 + self.expectedLoff = 2750.5 + self.expectedResolution = 2000. + self.expectedNblin = 5500 + self.expectedNbpix = 5500 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for in the filepath directory. + yaml_file = 'hima_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + # logger.info("Found myfiles = ", myfiles) + # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + # initHima() + + def initMtg(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_mtg_netcdficare(). + """ + self.netcdfName = tmp_path / "Mmultic1kmNC4_mtgi1_202406281000.nc" + self.filepath = tmp_path + + # self.buildNetcdf( + # self.netcdfName, 5568, "mtgi1", + # x0 = -5568000, dx = 2000.0000047, y0 = 5568000, cfac = 13642337, coff = 2784.5, + # listVisible= {"VIS004", "VIS005", "VIS006", "VIS008", "VIS009", + # "IR_013", "IR_016", "IR_022"}, + # listIR={"IR_038", "WV_062", "WV_073", "IR_087", + # "IR_097", "IR_105", "IR_123", "IR_133"}, + # coordinateSystemName = "GeosCoordinateSystem2km", + # nomImageNavigation = "ImageNavigation2km", + # nomX = "X2km", nomY = "Y2km", + # time_coverage_end = "2024-06-28T10:08:41Z365") + + listVisible = { + "VIS004", "VIS005", "VIS006", "VIS008", + "VIS009", "IR_013", "IR_016", "IR_022" + } + + self.buildNetcdf( + self.netcdfName, 11136, "mtgi1", + x0=-5568000, dx=1000.0000047, y0=5568000, + cfac=4.093316350596011E7, coff=5568.5, + listVisible=listVisible, listIR={"IR_038", "IR_105"}, + coordinateSystemName="GeosCoordinateSystem1km", + nomImageNavigation="ImageNavigation1km", + nomX="X1km", nomY="Y1km", time_coverage_end="2024-06-28T10:08:41Z365") + + # We will check that the parameters written in the dummy netcdf + # can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:08:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "Meteosat-12" + self.expectedSensor = "fci" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 4.093316350596011E7 # 13642337. + self.expectedLfac = 4.093316350596011E7 # 13642337. + self.expectedCoff = 5568.5 # 2784.5 + self.expectedLoff = 5568.5 # 2784.5 + self.expectedResolution = 1000. + self.expectedNblin = 11136 # 5568 + self.expectedNbpix = 11136 # 5568 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for in the filepath directory. + yaml_file = 'mtg_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + # logger.info("Found myfiles = ", myfiles) + # {'msg_netcdficare': ['/tmp/Mmultic3kmNC4_msg03_202406281000.nc']} + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + # initMtg() + + def initGoesr(self, tmp_path) : + """ + A fake netcdf is built. + A scene is built with the reader to be tested, applied to this netcdf. + Called by test_goesr_netcdficare(). + """ + self.netcdfName = tmp_path / "Emultic2kmNC4_goes16_202406281000.nc" + self.filepath = tmp_path + + listVisible = { + "VIS_004", "VIS_006", "VIS_008", "VIS_014", "VIS_016", "VIS_022" + } + listIR = { + "IR_039", "IR_062", "IR_069", "IR_073", "IR_085", + "IR_096", "IR_103", "IR_114", "IR_123", "IR_133" + } + self.buildNetcdf( + self.netcdfName, 5424, "goes16", + x0=-5434894.8, dx=2004.017288, y0=5434894.8, + cfac=20425338.9, coff=2712.5, + listVisible=listVisible, listIR=listIR, + coordinateSystemName="GeosCoordinateSystem2km", + nomImageNavigation="ImageNavigation2km", + nomX="X2km", nomY="Y2km", + time_coverage_end="2024-06-28T10:08:41.1Z") + + # We will check that the parameters written in the dummy netcdf + # can be read. + self.expectedStartTime = "2024-06-28T10:00:09" + self.expectedEndTime = "2024-06-28T10:08:41" + actualAltitude = 35786691 + 6378169 # 42164860.0 + actualLongitude = 0.1 + self.expectedPlatform = "GOES-16" + self.expectedSensor = "abi" + self.expectedAltitude = actualAltitude + self.expectedLongitude = actualLongitude + self.expectedCfac = 20425338.9 + self.expectedLfac = 20425338.9 + self.expectedCoff = 2712.5 + self.expectedLoff = 2712.5 + self.expectedResolution = 2000. + self.expectedNblin = 5424 + self.expectedNbpix = 5424 + + # To build a scene at the date 20240628_100000, + # a netcdf corresponding to msg_netcdficare + # is looked for in the filepath directory. + yaml_file = 'goesr_netcdficare' + myfiles = find_files_and_readers( + base_dir=self.filepath, start_time=datetime(2024, 6, 28, 10, 0), + end_time=datetime(2024, 6, 28, 10, 0), reader=yaml_file) + + self.scn = Scene(filenames=myfiles, reader=yaml_file) + # initGoesr() + + def buildNetcdf( + self, ncName, nbpix=3712, nomSatellite="msg03", + x0=-5570254, dx=3000.40604, y0=5570254, + cfac=1.3642337E7, coff=1857.0, + listVisible={}, listIR={}, + coordinateSystemName="GeosCoordinateSystem", + nomImageNavigation="ImageNavigation", + nomX="X", nomY="Y", + time_coverage_end="2024-06-28T10:12:41Z365") : """ ncName : tmp_path / Mmultic3kmNC4_msg03_202406281000.nc A dummy icare Meteo France netcdf is built here. - Called by init(). + Called by initMsg... + listVisible = {"VIS006", "VIS008", "IR_016"} + listIR = {"IR_039", "WV_062", "WV_073", "IR_087", "IR_097", "IR_108", + "IR_120", "IR_134"} """ if os.path.exists(ncName) : os.remove(ncName) @@ -201,11 +467,11 @@ def buildNetcdf(self, ncName) : ncName, mode="w", clobber=True, diskless=False, persist=False, format='NETCDF4') - ncfileOut.createDimension(u'ny', 3712) - ncfileOut.createDimension(u'nx', 3712) + ncfileOut.createDimension(u'ny', nbpix) + ncfileOut.createDimension(u'nx', nbpix) ncfileOut.createDimension(u'numerical_count', 65536) ncfileOut.setncattr("time_coverage_start", "2024-06-28T10:00:09Z383") - ncfileOut.setncattr("time_coverage_end", "2024-06-28T10:12:41Z365") + ncfileOut.setncattr("time_coverage_end", time_coverage_end) ncfileOut.setncattr("Area_of_acquisition", "globe") fill_value = -32768 @@ -214,7 +480,7 @@ def buildNetcdf(self, ncName) : shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr("id", "msg03") + var.setncattr("id", nomSatellite) var.setncattr("dst", 35786691.) var.setncattr("lon", float(0.1)) @@ -225,58 +491,61 @@ def buildNetcdf(self, ncName) : var.setncattr("longitude_of_projection_origin", 0.) var = ncfileOut.createVariable( - "GeosCoordinateSystem", "c", zlib=True, complevel=4, + coordinateSystemName, "c", zlib=True, complevel=4, shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr( - "GeoTransform", - "-5570254, 3000.40604, 0, 5570254, 0, -3000.40604") + + stringGeotransform = str(x0) + ", " + str(dx) + ", 0, " + stringGeotransform += str(-x0) + ", 0, " + str(-dx) + # -5570254, 3000.40604, 0, 5570254, 0, -3000.40604 + var.setncattr("GeoTransform", stringGeotransform) var = ncfileOut.createVariable( - "ImageNavigation", "c", zlib=True, complevel=4, + nomImageNavigation, "c", zlib=True, complevel=4, shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None) - var.setncattr("CFAC", 1.3642337E7) - var.setncattr("LFAC", 1.3642337E7) - var.setncattr("COFF", 1857.0) - var.setncattr("LOFF", 1857.0) + var.setncattr("CFAC", cfac) + var.setncattr("LFAC", cfac) + var.setncattr("COFF", coff) + var.setncattr("LOFF", coff) var = ncfileOut.createVariable( - "X", 'float32', u'nx', zlib=True, complevel=4, + nomX, 'float32', u'nx', zlib=True, complevel=4, shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None) - x0 = -5570254. - dx = 3000.40604 - var[:] = np.array(([(x0 + dx * i) for i in range(3712)])) + var[:] = np.array(([(x0 + dx * i) for i in range(nbpix)])) var = ncfileOut.createVariable( - "Y", 'float32', u'ny', zlib=True, complevel=4, + nomY, 'float32', u'ny', zlib=True, complevel=4, shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None) - y0 = 5570254. - dy = -3000.40604 - var[:] = np.array(([(y0 + dy * i) for i in range(3712)])) - - self.visibleChannelsCreation(ncfileOut, fill_value) - self.infrarougeChannelsCreation(ncfileOut, fill_value) - + y0 = -x0 + dy = -dx + var[:] = np.array(([(y0 + dy * i) for i in range(nbpix)])) + + self.visibleChannelsCreation( + ncfileOut, fill_value, listVisible, nbpix, coordinateSystemName) + self.infrarougeChannelsCreation( + ncfileOut, fill_value, listIR, nbpix, coordinateSystemName) ncfileOut.close # buildNetcdf() - def visibleChannelsCreation(self, ncfileOut, fill_value) : - for channel in {"VIS006", "VIS008", "IR_016"} : + def visibleChannelsCreation( + self, ncfileOut, fill_value, listVisible, nbpix, coordinateSystemName) : + + for channel in listVisible : var = ncfileOut.createVariable( channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None, fill_value=fill_value) var[:] = np.array( - ([[i * 2 for i in range(3712)] for j in range(3712)])) + ([[i * 2 for i in range(nbpix)] for j in range(nbpix)])) # Hundredths of albedo between 0 and 10000. var.setncattr("scale_factor", 0.01) var.setncattr("add_offset", 0.) var.setncattr("bandfactor", 20.76) - var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + var.setncattr("_CoordinateSystems", coordinateSystemName) var = ncfileOut.createVariable( "Albedo_to_Native_count_" + channel, 'short', @@ -287,24 +556,25 @@ def visibleChannelsCreation(self, ncfileOut, fill_value) : var[:] = np.array(([-9999 for i in range(65536)])) # In order to come back to the native datas on 10, 12 or 16 bits. - def infrarougeChannelsCreation(self, ncfileOut, fill_value) : - for channel in { - "IR_039", "WV_062", "WV_073", "IR_087", "IR_097", - "IR_108", "IR_120", "IR_134"} : + def infrarougeChannelsCreation( + self, ncfileOut, fill_value, listIR, nbpix, coordinateSystemName) : + + for channel in listIR : var = ncfileOut.createVariable( channel, 'short', ('ny', 'nx'), zlib=True, complevel=4, shuffle=True, fletcher32=False, contiguous=False, chunksizes=None, endian='native', least_significant_digit=None, fill_value=fill_value) var[:] = np.array( - ([[-9000 + j * 4 for i in range(3712)] for j in range(3712)])) + ([[-9000 + j * 4 for i in range(nbpix)] for j in range(nbpix)]) + ) # Hundredths of celcius degrees. var.setncattr("scale_factor", 0.01) var.setncattr("add_offset", 273.15) var.setncattr("nuc", 1600.548) var.setncattr("alpha", 0.9963) var.setncattr("beta", 2.185) - var.setncattr("_CoordinateSystems", "GeosCoordinateSystem") + var.setncattr("_CoordinateSystems", coordinateSystemName) var = ncfileOut.createVariable( "Temp_to_Native_count_" + channel, 'short', From df5452dfa488097822f33e519839dd35e640668c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?laurent=20p=C3=A9rier?= Date: Thu, 5 Dec 2024 18:51:00 +0100 Subject: [PATCH 9/9] Shorter test to understand what the matter is with the pull request. --- satpy/tests/reader_tests/test_geos_netcdficare.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/satpy/tests/reader_tests/test_geos_netcdficare.py b/satpy/tests/reader_tests/test_geos_netcdficare.py index 992971f44c..f0157e6ec3 100644 --- a/satpy/tests/reader_tests/test_geos_netcdficare.py +++ b/satpy/tests/reader_tests/test_geos_netcdficare.py @@ -52,7 +52,7 @@ class TestGeosNetcdfIcareReader() : # Test of the geos_netcdficare reader. # This reader has been build for the Icare Meteo France netcdfs. - + ''' def test_mtg_netcdficare(self, tmp_path) : """ A dummy netcdf is built. A scene self.scn for the nir_16 product for the same date @@ -77,7 +77,7 @@ def test_msg_netcdficare(self, tmp_path) : self.values = self.scn.values() self.checkingSceneParameter("convection") # test_msg_netcdficare() - + ''' def test_hima_netcdficare(self, tmp_path) : self.initHima(tmp_path) self.scn.load(['B10']) @@ -85,13 +85,14 @@ def test_hima_netcdficare(self, tmp_path) : self.checkingSceneParameter("B10") # test_hima_netcdficare() + ''' def test_goesr_netcdficare(self, tmp_path) : self.initGoesr(tmp_path) self.scn.load(['C02']) self.values = self.scn.values() self.checkingSceneParameter("C02") # test_goesr_netcdficare() - + ''' # ----------------------------------------------------- # # typeImage : convection, airmass... # ----------------------------------------------------- # @@ -356,7 +357,7 @@ def initMtg(self, tmp_path) : listVisible=listVisible, listIR={"IR_038", "IR_105"}, coordinateSystemName="GeosCoordinateSystem1km", nomImageNavigation="ImageNavigation1km", - nomX="X1km", nomY="Y1km", time_coverage_end="2024-06-28T10:08:41Z365") + nomX="X1km", nomY="Y1km", time_coverage_end="2024-06-28T10:08:41Z") # We will check that the parameters written in the dummy netcdf # can be read.