satpy icon indicating copy to clipboard operation
satpy copied to clipboard

Loss of UTM zone / CRS metadata when merging multiple Landsat-8 scenes

Open anikfal opened this issue 7 months ago • 3 comments

Single scene includes the correct area with a proj_dict (UTM zone 38N):

>>> type(scn_single)
<class 'xarray.core.dataarray.DataArray'>
>>> scn_single.attrs
{'AREA_OR_POINT': 'Point', '_FillValue': np.uint16(0), 'scale_factor': 1.0, 'add_offset': 0.0, 'perc_cloud_cover': 54.35, 'platform_name': 'Landsat-8', 'sensor': 'OLI_TIRS', 'standard_name': 'toa_bidirectional_reflectance', 'units': '%', 'saturated': True, 'start_time': datetime.datetime(2022, 5, 16, 7, 31, 57, tzinfo=datetime.timezone.utc), 'end_time': datetime.datetime(2022, 5, 16, 7, 31, 57, tzinfo=datetime.timezone.utc), 'reader': 'oli_tirs_l1_tif', 'area': Area ID: geotiff_area
Description: WGS84 / UTM zone 38N
Projection ID: WGS84 / UTM zone 38N
Projection: {'datum': 'WGS84', 'no_defs': 'None', 'proj': 'utm', 'type': 'crs', 'units': 'm', 'zone': '38'}
Number of columns: 7641
Number of rows: 7771
Area extent: (535185.0, 4190385.0, 764415.0, 4423515.0), 'name': 'B5', 'wavelength': WavelengthRange(min=0.845, central=0.867, max=0.885, unit='µm'), 'resolution': 30, 'calibration': 'reflectance', 'modifiers': (), '_satpy_id': DataID(name='B5', wavelength=WavelengthRange(min=0.845, central=0.867, max=0.885, unit='µm'), resolution=30, calibration=<1>, modifiers=()), 'ancillary_variables': []}

Merged multiple scenes retain UTM-like coordinates (x,y in meters), but no area or proj_dict is attached:

>>> type(scn_multiple)
<class 'xarray.core.dataarray.DataArray'>
>>> scn_multiple.attrs
{'add_offset': 0.0, 'AREA_OR_POINT': 'Point', '_FillValue': np.uint16(0), 'standard_name': 'toa_bidirectional_reflectance', 'units': '%', 'scale_factor': 1.0, 'sensor': 'OLI_TIRS', 'start_time': datetime.datetime(2022, 5, 16, 6, 42, 57, tzinfo=datetime.timezone.utc), 'end_time': datetime.datetime(2022, 5, 16, 7, 33, 9, tzinfo=datetime.timezone.utc), 'reader': 'oli_tirs_l1_tif', 'name': 'B5', 'wavelength': WavelengthRange(min=0.845, central=0.867, max=0.885, unit='µm'), 'resolution': 30, 'calibration': 'reflectance', 'modifiers': (), '_satpy_id': DataID(name='B5', wavelength=WavelengthRange(min=0.845, central=0.867, max=0.885, unit='µm'), resolution=30, calibration=<1>, modifiers=()), 'ancillary_variables': []}

This raises two concerns:

  1. If merged scenes span multiple UTM zones, the coordinates become inconsistent (the same easting/northing values have different meanings across zones).
  2. Even when all input scenes are from the same UTM zone, the loss of CRS metadata makes the merged coordinates ambiguous and unusable for resampling/reprojection without manual intervention.

anikfal avatar Aug 23 '25 17:08 anikfal

Is there a reason you didn't file a bug report and follow the template (https://github.com/pytroll/satpy/blob/main/.github%2FISSUE_TEMPLATE%2Fbug_report.md)? We'll need the code you're running to have any idea what's going on?

djhoese avatar Aug 23 '25 18:08 djhoese

And what version of satpy are you using?

djhoese avatar Aug 23 '25 18:08 djhoese

To be honest I'm not really sure what to do about this, if multiple UTM zones are crossed then I struggle to see an acceptable solution...

simonrp84 avatar Aug 25 '25 06:08 simonrp84