Loss of UTM zone / CRS metadata when merging multiple Landsat-8 scenes
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:
- If merged scenes span multiple UTM zones, the coordinates become inconsistent (the same easting/northing values have different meanings across zones).
- 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.
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?
And what version of satpy are you using?
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...