Skip to content
This repository was archived by the owner on Jun 2, 2025. It is now read-only.

Commit 92496f1

Browse files
authored
Add check for lat long location in satellite data spatial slicing (#350)
1 parent 9ec252e commit 92496f1

File tree

2 files changed

+50
-9
lines changed

2 files changed

+50
-9
lines changed

ocf_datapipes/select/select_spatial_slice.py

Lines changed: 17 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -103,24 +103,35 @@ def _get_idx_of_pixel_closest_to_poi(
103103

104104
def _get_idx_of_pixel_closest_to_poi_geostationary(
105105
xr_data: xr.DataArray,
106-
center_osgb: Location,
106+
center_coordinate: Location,
107107
) -> Location:
108108
"""
109109
Return x and y index location of pixel at center of region of interest.
110110
111111
Args:
112112
xr_data: Xarray dataset
113-
center_osgb: Center in OSGB coordinates
113+
center_coordinate: Central coordinate
114114
115115
Returns:
116116
Location for the center pixel in geostationary coordinates
117117
"""
118-
119118
xr_coords, xr_x_dim, xr_y_dim = spatial_coord_type(xr_data)
119+
if center_coordinate.coordinate_system == "osgb":
120+
x, y = osgb_to_geostationary_area_coords(
121+
x=center_coordinate.x, y=center_coordinate.y, xr_data=xr_data
122+
)
123+
elif center_coordinate.coordinate_system == "lon_lat":
124+
x, y = lon_lat_to_geostationary_area_coords(
125+
x=center_coordinate.x, y=center_coordinate.y, xr_data=xr_data
126+
)
127+
else:
128+
raise NotImplementedError(
129+
f"Only 'osgb' and 'lon_lat' location coordinates are \
130+
supported in conversion to geostationary \
131+
- not '{center_coordinate.coordinate_system}'"
132+
)
120133

121-
x, y = osgb_to_geostationary_area_coords(x=center_osgb.x, y=center_osgb.y, xr_data=xr_data)
122134
center_geostationary = Location(x=x, y=y, coordinate_system="geostationary")
123-
124135
# Check that the requested point lies within the data
125136
assert xr_data[xr_x_dim].min() < x < xr_data[xr_x_dim].max()
126137
assert xr_data[xr_y_dim].min() < y < xr_data[xr_y_dim].max()
@@ -390,7 +401,7 @@ def select_spatial_slice_pixels(
390401
if xr_coords == "geostationary":
391402
center_idx: Location = _get_idx_of_pixel_closest_to_poi_geostationary(
392403
xr_data=xr_data,
393-
center_osgb=location,
404+
center_coordinate=location,
394405
)
395406
else:
396407
center_idx: Location = _get_idx_of_pixel_closest_to_poi(

tests/select/test_select_spatial_slice.py

Lines changed: 33 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,16 @@
11
import numpy as np
22
import xarray as xr
3-
from ocf_datapipes.utils import Location
43

54
from ocf_datapipes.select import (
65
PickLocations,
76
SelectSpatialSliceMeters,
87
SelectSpatialSlicePixels,
98
)
10-
11-
from ocf_datapipes.select.select_spatial_slice import slice_spatial_pixel_window_from_xarray
9+
from ocf_datapipes.select.select_spatial_slice import (
10+
_get_idx_of_pixel_closest_to_poi_geostationary,
11+
slice_spatial_pixel_window_from_xarray,
12+
)
13+
from ocf_datapipes.utils import Location
1214

1315

1416
def test_slice_spatial_pixel_window_from_xarray_function():
@@ -158,3 +160,31 @@ def test_select_spatial_slice_meters_icon_global(passiv_datapipe, icon_global_da
158160
# ICON global has roughly 13km spacing, so this should be around 7x7 grid
159161
assert len(data.longitude) == 49
160162
assert len(data.latitude) == 49
163+
164+
165+
def test_get_idx_of_pixel_closest_to_poi_geostationary_lon_lat_location():
166+
# Create dummy data
167+
x = np.arange(5000000, -5000000, -5000)
168+
y = np.arange(5000000, -5000000, -5000)[::-1]
169+
170+
xr_data = xr.Dataset(
171+
data_vars=dict(
172+
data=(["x_geostationary", "y_geostationary"], np.random.normal(size=(len(x), len(y)))),
173+
),
174+
coords=dict(
175+
x_geostationary=(["x_geostationary"], x),
176+
y_geostationary=(["y_geostationary"], y),
177+
),
178+
)
179+
xr_data.attrs["area"] = (
180+
"msg_seviri_iodc_3km:\n description: MSG SEVIRI Indian Ocean Data Coverage service area definition with\n 3 km resolution\n projection:\n proj: geos\n lon_0: 41.5\n h: 35785831\n x_0: 0\n y_0: 0\n a: 6378169\n rf: 295.488065897014\n no_defs: null\n type: crs\n shape:\n height: 3712\n width: 3712\n area_extent:\n lower_left_xy: [5000000, 5000000]\n upper_right_xy: [-5000000, -5000000]\n units: m\n"
181+
)
182+
183+
center = Location(x=77.1, y=28.6, coordinate_system="lon_lat")
184+
185+
location_center_idx = _get_idx_of_pixel_closest_to_poi_geostationary(
186+
xr_data=xr_data, center_coordinate=center
187+
)
188+
189+
assert location_center_idx.coordinate_system == "idx"
190+
assert location_center_idx.x == 2000

0 commit comments

Comments
 (0)