Skip to content

Get the bounding box from the HDF EOS level-1 files and pass it to the array attributes #3399

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.

Sign up for GitHub

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

Draft
adybbroe wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Get the bounding box from the HDF EOS level-1 files and pass it to the array attributes #3399

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
055e299
Get the gringpoint (bounding box) data from the HDF EOS files, and pa…
May 21, 2026
a3e833b
Fix the unit tests
May 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion satpy/readers/core/hdfeos.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2019 Satpy developers
# Copyright (c) 2019-2026 Satpy developers
#
# This file is part of satpy.
#
Expand Down Expand Up @@ -122,6 +122,7 @@ def _load_all_metadata_attributes(self):
continue
else:
metadata.update(self.read_mda(str_val))

return metadata

@classmethod
Expand Down Expand Up @@ -225,6 +226,7 @@ def load_dataset(self, dataset_name, is_category=False):
dataset = self._read_dataset_in_file(dataset_name)
chunks = self._chunks_for_variable(dataset)
dask_arr = from_sds(dataset, self.filename, chunks=chunks)
#breakpoint()
dims = ("y", "x") if dask_arr.ndim == 2 else None
data = xr.DataArray(dask_arr, dims=dims,
attrs=dataset.attributes())
Expand Down Expand Up @@ -385,6 +387,11 @@ def get_dataset(self, dataset_id: DataID, dataset_info: dict) -> xr.DataArray:
data.attrs[key] = dataset_info[key]
self._add_satpy_metadata(dataset_id, data)

# mask = (~data.isnull().all(dim="x")).compute()
# idx = np.where(mask.values)[0]
# data_reduced = data.isel(y=idx)
# print(data_reduced.shape)

return data

def get_interpolated_dataset(self, dataset_name: str, resolution: int) -> xr.DataArray:
Expand Down
52 changes: 23 additions & 29 deletions satpy/readers/modis_l1b.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2010-2017 Satpy developers
# Copyright (c) 2010-2026 Satpy developers
#
# This file is part of satpy.
#
Expand Down Expand Up @@ -104,6 +104,14 @@ def __init__(self, filename, filename_info, filetype_info, mask_saturated=True,
super().__init__(filename, filename_info, filetype_info, **kwargs)
self._mask_saturated = mask_saturated

gringpoint = self.metadata["INVENTORYMETADATA"]["SPATIALDOMAINCONTAINER"][
"HORIZONTALSPATIALDOMAINCONTAINER"]["GPOLYGON"]["GPOLYGONCONTAINER"]["GRINGPOINT"]
bbox = {"longitudes": gringpoint["GRINGPOINTLONGITUDE"]["VALUE"],
"latitudes": gringpoint["GRINGPOINTLATITUDE"]["VALUE"]}
self._boundingbox = bbox
# Add the sub-satellite track geolocation as well if available? FIXME!
# Seems this is repeated for every HDF-EOS file being read.

ds = self.metadata["INVENTORYMETADATA"][
"COLLECTIONDESCRIPTIONCLASS"]["SHORTNAME"]["VALUE"]
self.resolution = self.res[ds[-3]]
Expand All @@ -126,37 +134,17 @@ def get_dataset(self, key, info):
array = self._fill_saturated(array, valid_max)
array = self._mask_invalid(array, valid_min, valid_max)
array = self._mask_uncertain_pixels(array, uncertainty, band_index)
# array = self._mask_rows_all_nan(array)
info["boundingbox"] = self._boundingbox
Copy link
Copy Markdown
Member

@djhoese djhoese May 21, 2026

Choose a reason for hiding this comment

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

"bbox"? And I wonder if this should go in the .attrs["orbital_parameters"]?

Copy link
Copy Markdown
Member

@djhoese djhoese May 21, 2026

Choose a reason for hiding this comment

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

Eh I guess it isn't "orbital" though. Nevermind. Ohand I just realized this is a draft. Sorry for early review.

Copy link
Copy Markdown
Contributor Author

@adybbroe adybbroe May 21, 2026

Choose a reason for hiding this comment

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

Do you prefer bbox over boundingbox? Or should it be bounding_box?

Copy link
Copy Markdown
Contributor Author

@adybbroe adybbroe May 21, 2026

Choose a reason for hiding this comment

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

...and yes "orbital_parameters" sounds reasonable

Copy link
Copy Markdown
Member

@djhoese djhoese May 21, 2026

Choose a reason for hiding this comment

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

I think I prefer bbox but bounding_box is maybe OK too. I think I take back what I said about orbital_parameters. I'm not sure it should go there. I think long term it will go in the .attrs of the future SwathDefinition, but the .attrs of the DataArray is good enough for now. But I also still think the min/max of all lon/lats is the way to fix this in pyresample. This bbox solution is a larger undertaking that needs discussion.

projectable = self._calibrate_data(key, info, array, var_attrs, band_index)

# if ((platform_name == 'Aqua' and key['name'] in ["6", "27", "36"]) or
# (platform_name == 'Terra' and key['name'] in ["29"])):
# height, width = projectable.shape
# row_indices = projectable.mask.sum(1) == width
# if row_indices.sum() != height:
# projectable.mask[row_indices, :] = True

# Get the orbit number
# if not satscene.orbit:
# mda = self.data.attributes()["CoreMetadata.0"]
# orbit_idx = mda.index("ORBITNUMBER")
# satscene.orbit = mda[orbit_idx + 111:orbit_idx + 116]

# Trimming out dead sensor lines (detectors) on terra:
# (in addition channel 27, 30, 34, 35, and 36 are nosiy)
# if satscene.satname == "terra":
# for band in ["29"]:
# if not satscene[band].is_loaded() or satscene[band].data.mask.all():
# continue
# width = satscene[band].data.shape[1]
# height = satscene[band].data.shape[0]
# indices = satscene[band].data.mask.sum(1) < width
# if indices.sum() == height:
# continue
# satscene[band] = satscene[band].data[indices, :]
# satscene[band].area = geometry.SwathDefinition(
# lons=satscene[band].area.lons[indices, :],
# lats=satscene[band].area.lats[indices, :])
self._add_satpy_metadata(key, projectable)

# mask = (~projectable.isnull().all(dim="x")).compute()
# idx = np.where(mask.values)[0]
# data_reduced = projectable.isel(y=idx)
# print(data_reduced.shape)

return projectable

def _get_band_variable_name_and_index(self, band_name):
Expand Down Expand Up @@ -213,11 +201,17 @@ def _mask_invalid(self, array, valid_min, valid_max):
def _mask_uncertain_pixels(self, array, uncertainty, band_index):
if not self._mask_saturated:
return array

uncertainty_chunks = self._chunks_for_variable(uncertainty)
band_uncertainty = from_sds(uncertainty, self.filename, chunks=uncertainty_chunks)[band_index, :, :]
array = array.where(band_uncertainty < 15)
return array

def _mask_rows_all_nan(self, array):
"""Mask rows where all values are NaN."""
arr_reduced = array.dropna(dim="y", how="all")
return arr_reduced

def _calibrate_data(self, key, info, array, var_attrs, index):
if key["calibration"] == "brightness_temperature":
projectable = calibrate_bt(array, var_attrs, index, key["name"])
Expand Down
18 changes: 16 additions & 2 deletions satpy/tests/reader_tests/modis_tests/_modis_fixtures.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2021 Satpy developers
# Copyright (c) 2021-2026 Satpy developers
#
# This file is part of satpy.
#
Expand Down Expand Up @@ -331,7 +331,21 @@ def _create_core_metadata(file_shortname: str) -> str:
f"VALUE = {file_shortname!r}\nEND_OBJECT = SHORTNAME\n\n" \
"OBJECT = VERSIONID\nNUM_VAL = 1\nVALUE = 6\nEND_OBJECT = VERSIONID\n\n" \
"END_GROUP = COLLECTIONDESCRIPTIONCLASS\n\n"
core_metadata_header += "\n\n" + inst_metadata + collection_metadata

gring_metadata = "GROUP = SPATIALDOMAINCONTAINER\n\n GROUP = HORIZONTALSPATIALDOMAINCONTAINER\n\n "\
'GROUP = GPOLYGON\n\n OBJECT = GPOLYGONCONTAINER\n CLASS = "1"\n\n '\
'GROUP = GRING\n CLASS = "1"\n\n OBJECT = EXCLUSIONGRINGFLAG\n NUM_VAL = 1\n '\
'CLASS = "1"\n VALUE = "N"\n END_OBJECT = EXCLUSIONGRINGFLAG\n\n END_GROUP = GRING\n\n '\
'GROUP = GRINGPOINT\n CLASS = "1"\n\n OBJECT = GRINGPOINTLONGITUDE\n NUM_VAL = 4\n '\
'CLASS = "1"\n VALUE = (37.5053242268888, 6.38401236692491, -48.2870137383305, 36.7762332873298)\n '\
'END_OBJECT = GRINGPOINTLONGITUDE\n\n OBJECT = GRINGPOINTLATITUDE\n NUM_VAL = 4\n CLASS = "1"\n '\
"VALUE = (34.5168084728086, 30.695318900416, 70.1120029478276, 83.8447003025456)\n "\
"END_OBJECT = GRINGPOINTLATITUDE\n\n OBJECT = GRINGPOINTSEQUENCENO\n NUM_VAL = 4\n "\
'CLASS = "1"\n VALUE = (1, 2, 3, 4)\n END_OBJECT = GRINGPOINTSEQUENCENO\n\n '\
"END_GROUP = GRINGPOINT\n\n END_OBJECT = GPOLYGONCONTAINER\n\n END_GROUP = GPOLYGON\n\n "\
"END_GROUP = HORIZONTALSPATIALDOMAINCONTAINER\n\n END_GROUP = SPATIALDOMAINCONTAINER\n\n"

core_metadata_header += "\n\n" + inst_metadata + collection_metadata + gring_metadata
return core_metadata_header


Expand Down
Loading