Skip to content
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

Jump to bottom

Fix geostationary bbox having inf values #493

Merged
merged 7 commits into from
Feb 7, 2023
Merged

Fix geostationary bbox having inf values #493

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
9a9a16f
Fix geostationary bbox having inf values
mraspaud Jan 26, 2023
dd2ed80
Fix geostationary bbox having inf values
mraspaud Feb 3, 2023
493e859
Fix goes bbox computation for gradient search tests
mraspaud Feb 3, 2023
6b04f6d
Remove unneeded comments
mraspaud Feb 3, 2023
33553c5
Add shapely as a hard dependency
mraspaud Feb 6, 2023
093e706
Split utils import
mraspaud Feb 6, 2023
b20f66c
Remove 32-bit deploy
mraspaud Feb 7, 2023
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
49 changes: 29 additions & 20 deletions pyresample/geometry.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,11 +38,11 @@
from pyresample.boundary import AreaDefBoundary, Boundary, SimpleBoundary
from pyresample.utils import (
check_slice_orientation,
get_geostationary_height,
load_cf_area,
proj4_dict_to_str,
proj4_radius_parameters,
)
from pyresample.utils.proj4 import get_geostationary_height
Copy link
Member

Choose a reason for hiding this comment

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

Just curious, why the import change?

Copy link
Member Author

Choose a reason for hiding this comment

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

I think it's more explicit this way

Copy link
Member

Choose a reason for hiding this comment

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

Sure, but maybe less consistent? Maybe all the utils should be explicit and we gret of the convenience import in utils/__init__.py?


try:
from xarray import DataArray
Expand Down Expand Up @@ -1556,7 +1556,7 @@ def _get_geo_boundary_sides(self, frequency=None):
# Ensure an even number of vertices for side creation
if (frequency % 2) != 0:
frequency = frequency + 1
lons, lats = get_geostationary_bounding_box(self, nb_points=frequency)
lons, lats = get_geostationary_bounding_box_in_lonlats(self, nb_points=frequency)
# Retrieve dummy sides for GEO (side1 and side3 always of length 2)
side02_step = int(frequency / 2) - 1
lon_sides = [lons[slice(0, side02_step + 1)],
Expand Down Expand Up @@ -2772,40 +2772,49 @@ def get_geostationary_angle_extent(geos_area):
aeq = 1 - req ** 2 / (h ** 2)
ap_ = 1 - rp ** 2 / (h ** 2)

# generate points around the north hemisphere in satellite projection
# make it a bit smaller so that we stay inside the valid area
xmax = np.arccos(np.sqrt(aeq))
ymax = np.arccos(np.sqrt(ap_))
return xmax, ymax
x_angle = np.arccos(np.sqrt(aeq))
y_angle = np.arccos(np.sqrt(ap_))
return x_angle, y_angle


def get_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
"""Get the bbox in geos projection coordinates of the valid pixels inside `geos_area`.

Notes:
- The first and last element of the output vectors are equal.
- If nb_points is even, it will return x and y vectors of length nb_points + 1.
Comment on lines -2785 to -2787
Copy link
Member

Choose a reason for hiding this comment

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

Is it a problem for any of the other interfaces (the bounds and polygon stuff) that the end points are not the same?

Copy link
Member Author

Choose a reason for hiding this comment

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

The tests didn't break, so I think we're good. Also I always thought that was the way we did it in pyresample.


Parameters
----------
nb_points : Number of points on the polygon.

"""
xmax, ymax = get_geostationary_angle_extent(geos_area)
x, y = get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points)
ll_x, ll_y, ur_x, ur_y = geos_area.area_extent

from shapely.geometry import Polygon
geo_bbox = Polygon(np.vstack((x, y)).T)
area_bbox = Polygon(((ll_x, ll_y), (ll_x, ur_y), (ur_x, ur_y), (ur_x, ll_y)))
intersection = area_bbox.intersection(geo_bbox)
Comment on lines +2790 to +2793
Copy link
Member

Choose a reason for hiding this comment

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

There is a very good chance I'm missing something, but could this be simplified (logically, maybe not in code) to min/max each of the coordinates in the full geostationary bounding box with the extents of the area? So if the full bbox says a point is at (10_000, -85_000_000) and the extents of the area are (-10_000, -60_000_000, 50_000, 10_000_000) then we know that the full bbox point's x coordinate is within the extent, but the y coordinate can be clipped to -60_000_000.

Not sure this is accurate, or makes sense, or is a good idea, but it is the first thing that comes to mind.

Copy link
Member Author

Choose a reason for hiding this comment

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

That's basically what we had before, in some way. The problem is to get the corner points of the intersection between the disc and the truncation, especially if only one corner of the truncated area is in space.

Copy link
Member

Choose a reason for hiding this comment

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

Ok so you're saying if the disk polygon was computed and has points partially in space that clipping to the extents of the area wouldn't effect some of those points? How does this current intersection do that?

I thought the problem in main was that the calculation of the full disk polygon coordinates were including points outside the full disk? I'll try reading the code to see if I can understand it more.

In your comment above you said:

The problem appeared when we just had an area describing a portion of the earth that did not cover the equator. In that case, the simple clipping we did to "cut" the disk was giving (unneeded) points outside the valid disk.

When you say "cut", do you mean clipping the disk coordinates to the extent of the area or do you mean how the coordinates were calculated for half the disk and then replicated for the other half?

Copy link
Member Author

Choose a reason for hiding this comment

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

ok, so the polygon for the full disk looks like this:
full_disk

Now, if we just clip y at 2e6 like it is in main, we get this:
clipped_y

If you look carefully, you see the last point on each side of the y=2e6 line ends up in space, because we 'project' the extreme points of the equator onto a latitude where the sphere is not as wide.

So how do we avoid this?
One solution would be to not compute everything below y=2e6 and find the edge of the disk at that y position to close the polygon. It's not too hard, but there are many different case, whether the disk is cut horizontally, vertically, or both, includes the equator or prime meridian, or not.
A simpler geometrical approach in my opinion is just to find the intersection of the full disk polygon with the rectangle of the area extent. Since they are all in the same projection, we can see that as a simple 2d geometry problem.
full_and_box

Copy link
Member

@djhoese djhoese Feb 3, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

Ok, I think I got it. One more alternative solution: mask out (remove) any points not within the bbox. The main point here is that our extent bbox is not a complicated polygon so we only have two limits to check:

disk_points = ...
area_points = disk_points[
    (disk_points[:, 0] >= min(extent[0], extent[2])) &
    (disk_points[:, 0] <= max(extent[0], extent[2])) &
    (disk_points[:, 1] >= min(extent[1], extent[3])) &
    (disk_points[:, 1] <= min(extent[1], extent[3]))]

Copy link
Member

Choose a reason for hiding this comment

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

The only other thing, which is at this point getting ridiculous and we should just use shapely, would be to combine my two suggestions and every outer disk pixel gets clipped to the min/max of the inner disk pixels (the pixels inside the area's bbox).

Copy link
Member

Choose a reason for hiding this comment

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

Nevermind again, this last suggestion doesn't fix the issue. You still end up with the points being no closer to the extent of the area. Ok, shapely it is. Maybe at this point we say shapely is required? What is it being used for now? Did @ghiggi want to depend on it more heavily in the future for spherical stuff?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, I plan to use it for several tasks. Could make sense to have shapely as a requirement.
Just check if we need to constrain to shapely < 2.0.0 or we ask for shapely > 2.0.0

Copy link
Member

Choose a reason for hiding this comment

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

Maybe at this point we say shapely is required? What is it being used for now?

Gradient search resampler is using shapely to find the overlapping source and target area chunk pairs so that only those are handled.

Copy link
Member Author

Choose a reason for hiding this comment

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

Made shapely a hard requirement

try:
x, y = intersection.boundary.xy
except NotImplementedError:
return [], []
return np.asanyarray(x[:-1]), np.asanyarray(y[:-1])


def get_full_geostationary_bounding_box_in_proj_coords(geos_area, nb_points=50):
"""Get the bbox in geos projection coordinates of the full disk in `geos_area` projection.

Args:
nb_points: Number of points on the polygon
"""
x_max_angle, y_max_angle = get_geostationary_angle_extent(geos_area)
h = get_geostationary_height(geos_area.crs)

# generate points around the north hemisphere in satellite projection
# make it a bit smaller so that we stay inside the valid area
x = np.cos(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (xmax - 0.0001)
y = -np.sin(np.linspace(-np.pi, 0, int(nb_points / 2.0) + 1)) * (ymax - 0.0001)

ll_x, ll_y, ur_x, ur_y = geos_area.area_extent

points_around = np.linspace(-np.pi, np.pi, nb_points, endpoint=False)
x = np.cos(points_around) * (x_max_angle - 0.0001)
y = -np.sin(points_around) * (y_max_angle - 0.0001)
x *= h
y *= h
# We remove one element with [:-1] to avoid duplicate values at the equator
x = np.clip(np.concatenate([x[:-1], x[::-1]]), min(ll_x, ur_x), max(ll_x, ur_x))
y = np.clip(np.concatenate([y[:-1], -y]), min(ll_y, ur_y), max(ll_y, ur_y))

return x, y

Expand Down
4 changes: 2 additions & 2 deletions pyresample/gradient/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
from pyresample.geometry import (
AreaDefinition,
SwathDefinition,
get_geostationary_bounding_box,
get_geostationary_bounding_box_in_lonlats,
)
from pyresample.gradient._gradient_search import (
one_step_gradient_indices,
Expand Down Expand Up @@ -375,7 +375,7 @@ def _check_input_coordinates(dst_x, dst_y,
def get_border_lonlats(geo_def):
"""Get the border x- and y-coordinates."""
if geo_def.proj_dict['proj'] == 'geos':
lon_b, lat_b = get_geostationary_bounding_box(geo_def, 3600)
lon_b, lat_b = get_geostationary_bounding_box_in_lonlats(geo_def, 3600)
else:
lons, lats = geo_def.get_boundary_lonlats()
lon_b = np.concatenate((lons.side1, lons.side2, lons.side3, lons.side4))
Expand Down
146 changes: 119 additions & 27 deletions pyresample/test/test_geometry.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2591,8 +2591,100 @@ def _get_fake_antimeridian_lonlats(use_dask: bool) -> tuple:
return lons, lats


class TestCrop(unittest.TestCase):
"""Test the area helpers."""
@pytest.fixture
def truncated_geos_area():
"""Create a truncated geostationary area."""
projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos',
'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
area_extent = (5567248.0742, 5570248.4773, -5570248.4773, 1393687.2705)
width = 3712
height = 1392
geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent)
return geos_area


@pytest.fixture
def truncated_geos_area_in_space():
"""Create a truncated geostationary area."""
projection = {'a': '6378169', 'h': '35785831', 'lon_0': '9.5', 'no_defs': 'None', 'proj': 'geos',
'rf': '295.488065897014', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
area_extent = (5575000, 5575000, 5570000, 5570000)
width = 10
height = 10
geos_area = geometry.AreaDefinition('msg_rss', "msg_rss", "msg_rss", projection, width, height, area_extent)
return geos_area


class TestGeostationaryTools:
"""Test the geostationary bbox tools."""

def test_get_full_geostationary_bbox(self, truncated_geos_area):
nb_points = 20
x, y = geometry.get_full_geostationary_bounding_box_in_proj_coords(truncated_geos_area, nb_points)
assert len(x) == nb_points
assert len(y) == nb_points

assert x[0] != x[-1]
assert y[0] != y[-1]

expected_x = np.array([-5.43062255e+06, -5.16482897e+06, -4.39346593e+06, -3.19203985e+06,
-1.67815466e+06, 3.32529726e-10, 1.67815466e+06, 3.19203985e+06,
4.39346593e+06, 5.16482897e+06, 5.43062255e+06, 5.16482897e+06,
4.39346593e+06, 3.19203985e+06, 1.67815466e+06, 3.32529726e-10,
-1.67815466e+06, -3.19203985e+06, -4.39346593e+06, -5.16482897e+06])

expected_y = np.array([6.62789871e-10, 1.67242779e+06, 3.18114670e+06, 4.37847280e+06,
5.14720348e+06, 5.41209002e+06, 5.14720348e+06, 4.37847280e+06,
3.18114670e+06, 1.67242779e+06, -0.00000000e+00, -1.67242779e+06,
-3.18114670e+06, -4.37847280e+06, -5.14720348e+06, -5.41209002e+06,
-5.14720348e+06, -4.37847280e+06, -3.18114670e+06, -1.67242779e+06])

np.testing.assert_allclose(x, expected_x)
np.testing.assert_allclose(y, expected_y)

def test_get_geostationary_bbox_works_with_truncated_area(self, truncated_geos_area):
"""Ensure the geostationary bbox works when truncated."""
lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)

expected_lon = np.array(
[-64.24072434653284, -68.69662326361153, -65.92516214783112, -60.726360278290336,
-47.39851775032484, 9.500000000000018, 66.39851775032487, 79.72636027829033,
84.92516214783113, 87.69662326361151, 83.24072434653286])
expected_lat = np.array(
[14.554922655532085, 17.768795771961937, 35.34328897185421, 52.597860701318254, 69.00533141646078,
79.1481121862375, 69.00533141646076, 52.597860701318254, 35.34328897185421, 17.768795771961933,
14.554922655532085])
np.testing.assert_allclose(lon, expected_lon)
np.testing.assert_allclose(lat, expected_lat)

def test_get_geostationary_bbox_works_with_truncated_area_proj_coords(self, truncated_geos_area):
"""Ensure the geostationary bbox works when truncated."""
x, y = geometry.get_geostationary_bounding_box_in_proj_coords(truncated_geos_area, 20)

expected_x = np.array(
[-5209128.302753595, -5164828.965702432, -4393465.934674804, -3192039.8468840676, -1678154.6586309497,
3.325297262895822e-10, 1678154.6586309501, 3192039.846884068, 4393465.934674805, 5164828.965702432,
5209128.302753594])
expected_y = np.array(
[1393687.2705, 1672427.7900638399, 3181146.6955466354, 4378472.798117005, 5147203.47659387,
5412090.016106332, 5147203.476593869, 4378472.798117005, 3181146.695546635, 1672427.7900638392,
1393687.2705])

np.testing.assert_allclose(x, expected_x)
np.testing.assert_allclose(y, expected_y)

def test_get_geostationary_bbox_does_not_contain_inf(self, truncated_geos_area):
"""Ensure the geostationary bbox does not contain np.inf."""
lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area, 20)
assert not any(np.isinf(lon))
assert not any(np.isinf(lat))

def test_get_geostationary_bbox_returns_empty_lonlats_in_space(self, truncated_geos_area_in_space):
"""Ensure the geostationary bbox is empty when in space."""
lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(truncated_geos_area_in_space, 20)

assert len(lon) == 0
assert len(lat) == 0

def test_get_geostationary_bbox(self):
"""Get the geostationary bbox."""
Expand All @@ -2606,23 +2698,20 @@ def test_get_geostationary_bbox(self):
geos_area.crs = CRS(proj_dict)
geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]

lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
# This musk be equal to lon.
elon = np.array([-79.23372832, -78.19662326, -75.42516215, -70.22636028,
-56.89851775, 0., 56.89851775, 70.22636028,
75.42516215, 78.19662326, 79.23372832, 78.19662326,
75.42516215, 70.22636028, 56.89851775, 0.,
-56.89851775, -70.22636028, -75.42516215, -78.19662326,
-79.23372832])
elat = np.array([0., 17.76879577, 35.34328897, 52.5978607,
69.00533142, 79.14811219, 69.00533142, 52.5978607,
35.34328897, 17.76879577, -0., -17.76879577,
-35.34328897, -52.5978607, -69.00533142, -79.14811219,
-69.00533142, -52.5978607, -35.34328897, -17.76879577,
0.])

np.testing.assert_allclose(lon, elon, atol=1e-07)
np.testing.assert_allclose(lat, elat, atol=1e-07)
lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
expected_lon = np.array([-78.19662326, -75.42516215, -70.22636028,
-56.89851775, 0., 56.89851775, 70.22636028,
75.42516215, 78.19662326, 79.23372832, 78.19662326,
75.42516215, 70.22636028, 56.89851775, 0.,
-56.89851775, -70.22636028, -75.42516215, -78.19662326, -79.23372832, ])
expected_lat = np.array([17.76879577, 35.34328897, 52.5978607,
69.00533142, 79.14811219, 69.00533142, 52.5978607,
35.34328897, 17.76879577, -0., -17.76879577,
-35.34328897, -52.5978607, -69.00533142, -79.14811219,
-69.00533142, -52.5978607, -35.34328897, -17.76879577, 0.])

np.testing.assert_allclose(lon, expected_lon, atol=1e-07)
np.testing.assert_allclose(lat, expected_lat, atol=1e-07)

geos_area = MagicMock()
lon_0 = 10
Expand All @@ -2634,8 +2723,8 @@ def test_get_geostationary_bbox(self):
geos_area.crs = CRS(proj_dict)
geos_area.area_extent = [-5500000., -5500000., 5500000., 5500000.]

lon, lat = geometry.get_geostationary_bounding_box(geos_area, 20)
np.testing.assert_allclose(lon, elon + lon_0)
lon, lat = geometry.get_geostationary_bounding_box_in_lonlats(geos_area, 20)
np.testing.assert_allclose(lon, expected_lon + lon_0)

def test_get_geostationary_angle_extent(self):
"""Get max geostationary angles."""
Expand Down Expand Up @@ -2675,6 +2764,10 @@ def test_get_geostationary_angle_extent(self):
np.testing.assert_allclose(expected,
geometry.get_geostationary_angle_extent(geos_area))


class TestCrop(unittest.TestCase):
"""Test the area helpers."""

def test_sub_area(self):
"""Sub area slicing."""
area = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
Expand Down Expand Up @@ -3156,20 +3249,19 @@ def test_geostationary_projection(self):
# Check boundary vertices
n_vertices = 10
boundary = areadef.boundary(frequency=n_vertices, force_clockwise=False)
boundary.vertices.shape

# Check boundary vertices is in correct order
expected_vertices = np.array([[-7.92337283e+01, 6.94302533e-15],
[-7.54251621e+01, 3.53432890e+01],
expected_vertices = np.array([[-7.54251621e+01, 3.53432890e+01],
[-5.68985178e+01, 6.90053314e+01],
[5.68985178e+01, 6.90053314e+01],
[7.54251621e+01, 3.53432890e+01],
[7.92337283e+01, -6.94302533e-15],
[7.92337283e+01, -0.00000000e+00],
[7.54251621e+01, -3.53432890e+01],
[5.68985178e+01, -6.90053314e+01],
[-5.68985178e+01, -6.90053314e+01],
[-7.54251621e+01, -3.53432890e+01]])
assert np.allclose(expected_vertices, boundary.vertices)
[-7.54251621e+01, -3.53432890e+01],
[-7.92337283e+01, 6.94302533e-15]])
np.testing.assert_allclose(expected_vertices, boundary.vertices)

def test_global_platee_caree_projection(self):
"""Test boundary for global platee caree projection."""
Expand Down
2 changes: 1 addition & 1 deletion pyresample/test/test_gradient.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -599,7 +599,7 @@ def test_check_overlap():
assert check_overlap(poly1, poly2) is False


@mock.patch('pyresample.gradient.get_geostationary_bounding_box')
@mock.patch('pyresample.gradient.get_geostationary_bounding_box_in_lonlats')
Copy link
Member

Choose a reason for hiding this comment

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

I really don't like how all these old tests mock the functions being called. I know this isn't what this PR is about, but mocking the bounding box function called by the border function in this test, what's being tested then? The same for the test below it, lots of mocks. In the future it'd be nice if we could construct an area definition that just produced the results we wanted.

Copy link
Member Author

Choose a reason for hiding this comment

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

I agree. Let's aim to keep mocking at a minimum.

def test_get_border_lonlats(get_geostationary_bounding_box):
"""Test that correct methods are called in get_border_lonlats()."""
from pyresample.gradient import get_border_lonlats
Expand Down
9 changes: 4 additions & 5 deletions setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,20 +34,19 @@
if sys.version_info < (3, 10):
requirements.append('importlib_metadata')

test_requires = ['rasterio', 'dask', 'xarray', 'cartopy', 'pillow', 'matplotlib', 'scipy', 'zarr',
'pytest-lazy-fixtures']
test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr',
'pytest-lazy-fixtures', 'shapely']
extras_require = {'numexpr': ['numexpr'],
'quicklook': ['matplotlib', 'cartopy>=0.20.0', 'pillow'],
'rasterio': ['rasterio'],
'dask': ['dask>=0.16.1'],
'cf': ['xarray'],
'gradient_search': ['shapely'],
'xarray_bilinear': ['xarray', 'dask', 'zarr'],
'tests': test_requires}
'tests': test_requires,
"geos_areas": ["shapely"]}

setup_requires = ['numpy>=1.10.0', 'cython']
test_requires = ['rasterio', 'dask', 'xarray', 'cartopy>=0.20.0', 'pillow', 'matplotlib', 'scipy', 'zarr',
'pytest-lazy-fixture']

if sys.platform.startswith("win"):
extra_compile_args = []
Expand Down