Skip to content

Add reader for EarthCARE MSI L1 data #3080

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open

Conversation

pdebuyl
Copy link
Contributor

@pdebuyl pdebuyl commented Mar 14, 2025

This PR adds a reader to satpy for EarthCARE MSI L1 data. It was initiated by @simonrp84

  • Tests added
  • Fully documented

Copy link

codecov bot commented Mar 14, 2025

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.19%. Comparing base (f15bd3b) to head (3f41f5b).
Report is 19 commits behind head on main.

Additional details and impacted files
@@           Coverage Diff            @@
##             main    #3080    +/-   ##
========================================
  Coverage   96.18%   96.19%            
========================================
  Files         391      393     +2     
  Lines       56633    56734   +101     
========================================
+ Hits        54471    54573   +102     
+ Misses       2162     2161     -1     
Flag Coverage Δ
behaviourtests 3.82% <0.00%> (-0.01%) ⬇️
unittests 96.28% <100.00%> (+<0.01%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@coveralls
Copy link

coveralls commented Mar 14, 2025

Pull Request Test Coverage Report for Build 14446587353

Details

  • 102 of 102 (100.0%) changed or added relevant lines in 2 files are covered.
  • No unchanged relevant lines lost coverage.
  • Overall coverage increased (+0.008%) to 96.294%

Totals Coverage Status
Change from base Build 14401317584: 0.008%
Covered Lines: 54828
Relevant Lines: 56938

💛 - Coveralls

@pdebuyl
Copy link
Contributor Author

pdebuyl commented Mar 18, 2025

Here's a self-contained example if it helps in reviewing :-)

Call as python show_one_MSI.py ~/nobackup/ECA/ECA_EXAF_MSI_RGR_1C_20250318T092816Z_20250318T100351Z_04566A.h5 --band TIR3

When given several input files, they are shown together.

import argparse
import matplotlib.pyplot as plt
import satpy
import cartopy
from satpy.writers import get_enhanced_image
import numpy as np

ECA_MSI_BANDS = ['SWIR1', 'SWIR2', 'TIR1', 'TIR2', 'TIR3', 'VIS', 'VNIR']

parser = argparse.ArgumentParser()
parser.add_argument('input', nargs='+')
parser.add_argument('--band', required=True, choices=ECA_MSI_BANDS)
parser.add_argument('-v', action="store_true")
args = parser.parse_args()


fig = plt.figure(figsize=(5,3.25))
crs = cartopy.crs.PlateCarree()
ax = fig.add_subplot(111, projection=crs)
ax.coastlines()
ax.gridlines(draw_labels=True)

scn = satpy.Scene(filenames=args.input, reader='msi_l1c_earthcare')
scn.load([args.band, 'latitude', 'longitude'])

if args.band in scn.available_dataset_names():
    # For bands, simply load calibrated data
    im = scn[args.band][:,:]
else:
    # For composites, use get_enhanced_image and convert to RGB
    im = get_enhanced_image(scn[args.band])
    im.stretch()
    tmp = (im.data.transpose('y', 'x', 'bands')*255).astype(np.uint8)
    im = tmp

lats = scn['latitude'].data
lons = scn['longitude'].data
ax.pcolormesh(lons, lats, im, transform=crs)

plt.show()

Comment on lines 22 to 84
class FakeHDF5FileHandler2(FakeHDF5FileHandler):
"""Swap-in HDF5 File Handler."""


def _setup_test_data(self, N_BANDS, N_SCANS, N_COLS):
# Set some default attributes
data = {
"ScienceData/pixel_values":
xr.DataArray(
da.ones((N_BANDS, N_SCANS, N_COLS), chunks=1024, dtype=np.float32),
attrs={"units": "Wm-2 sr-1 or K", "DIMENSION_LIST": DIMLIST},
dims=("band", "dim_2", "dim_1")),
"ScienceData/land_flag":
xr.DataArray(
da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.uint16),
attrs={"units": ""},
dims=("along_track", "across_track")),
"ScienceData/solar_azimuth_angle":
xr.DataArray(
da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.float32),
attrs={"units": "degrees"},
dims=("along_track", "across_track")),
"ScienceData/longitude":
xr.DataArray(
da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.float32),
attrs={"units": "degrees"},
dims=("along_track", "across_track")),
"ScienceData/latitude":
xr.DataArray(
da.ones((N_SCANS, N_COLS), chunks=1024, dtype=np.float32),
attrs={"units": "degrees"},
dims=("along_track", "across_track")),
"ScienceData/solar_spectral_irradiance":
xr.DataArray(
da.array(SOL_IRRAD),
attrs={"units": "W m-2"},
dims=("band", "across_track")),
}

return data

def get_test_content(self, filename, filename_info, filetype_info):
"""Mimic reader input file content."""
test_content = self._setup_test_data(N_BANDS, N_SCANS, N_COLS)
return test_content


class ECMSIL1CTester:
"""Test MSI/EarthCARE L1C Reader."""

def setup_method(self):
"""Wrap HDF5 file handler with our own fake handler."""
from satpy._config import config_search_paths
from satpy.readers.msi_ec_l1c_h5 import MSIECL1CFileHandler
self.reader_configs = config_search_paths(os.path.join("readers", self.yaml_file))
# http://stackoverflow.com/questions/12219967/how-to-mock-a-base-class-with-python-mock-library
self.p = mock.patch.object(MSIECL1CFileHandler, "__bases__", (FakeHDF5FileHandler2,))
self.fake_handler = self.p.start()
self.p.is_local = True

def teardown_method(self):
"""Stop wrapping the HDF5 file handler."""
self.p.stop()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nowadays we try to create stub files using fixtures instead of mocking, to you think this is something we could do here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you point me to a test doing that? It is often confusing to know which way is "the good way" :-)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here an example, I hope it's good enough

@pytest.fixture(scope="session")
def insat_filename(tmp_path_factory):
"""Create a fake insat 3d l1b file."""
filename = tmp_path_factory.mktemp("data") / "3DIMG_25OCT2022_0400_L1B_STD_V01R00.h5"
with h5netcdf.File(filename, mode="w") as h5f:
h5f.dimensions = dimensions
h5f.attrs.update(global_attrs)
for resolution, channels in CHANNELS_BY_RESOLUTION.items():
_create_channels(channels, h5f, resolution)
_create_lonlats(h5f, resolution)
return filename

:D

@mraspaud
Copy link
Member

Could you paste the image of a composite here? It’s always nice to have as a reference.

@simonrp84
Copy link
Member

If Pierre is busy, I'll add an example image this evening :)

@pdebuyl
Copy link
Contributor Author

pdebuyl commented Apr 6, 2025

One VIS (orbit file 04864D)

ECA_MSI_04864D_VIS

and ash (same file)

ECA_MSI_04864D_ash

Note: when plotting the ash composite, I often have this black region filling the space between the data and the side of the image.

@pdebuyl
Copy link
Contributor Author

pdebuyl commented Apr 14, 2025

Hello @simonrp84

The tests run fine but the decreasing coverage makes me realize that I don't use the DIMENSIONLIST attribute. So, I never user _fix_dims. I don't know whether it is necessary though, as data = data[band_index] should match the proper dimension.

@pdebuyl
Copy link
Contributor Author

pdebuyl commented Apr 14, 2025

Note (for Simon): I didn't actually change the reader itself, so all the data loads "as before". Also: the product user guide does not mention any DIMENSION_LIST attribute -> how does it appear ?

@simonrp84
Copy link
Member

I suspect that may have been some artefact in the commissioning data that was removed.
If everything loads nicely without it, then no problem and we don't need that any more!

The code did not seem necessary for production data.
@pdebuyl
Copy link
Contributor Author

pdebuyl commented Apr 14, 2025

Hi Simon, thanks for checking. I removed the code. I also checked that I could actually load MSI data from disk :-)

@pdebuyl
Copy link
Contributor Author

pdebuyl commented Apr 14, 2025

ok @mraspaud the tests are now file-based. only ubuntu latest fails but I think that this is the case for other open PRs.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants