Skip to content

Commit

Permalink
Merge pull request #342 from djhoese/bugfix-dynamic-area-def-antimeri…
Browse files Browse the repository at this point in the history
…dian

Fix DynamicAreaDefinition not handling lons over antimeridian
  • Loading branch information
djhoese authored Mar 22, 2021
2 parents 68414d7 + 5e3861a commit 9801978
Show file tree
Hide file tree
Showing 2 changed files with 116 additions and 59 deletions.
103 changes: 65 additions & 38 deletions pyresample/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -816,34 +816,43 @@ class DynamicAreaDefinition(object):
of the area to a given set of longitudes and latitudes, such that e.g.
polar satellite granules can be resampled optimally to a given projection.
Parameters
----------
area_id:
The name of the area.
description:
The description of the area.
projection:
The dictionary or string or CRS object of projection parameters.
Doesn't have to be complete. If not complete, ``proj_info`` must
be provided to ``freeze`` to "fill in" any missing parameters.
width:
x dimension in number of pixels, aka number of grid columns
height:
y dimension in number of pixels, aka number of grid rows
shape:
Corresponding array shape as (height, width)
area_extent:
The area extent of the area.
pixel_size_x:
Pixel width in projection units
pixel_size_y:
Pixel height in projection units
resolution:
Resolution of the resulting area as (pixel_size_x, pixel_size_y) or a scalar if pixel_size_x == pixel_size_y.
optimize_projection:
Whether the projection parameters have to be optimized.
rotation:
Rotation in degrees (negative is cw)
Note that if the provided projection is geographic (lon/lat degrees) and
the provided longitude and latitude data crosses the anti-meridian
(-180/180), the resulting area will be the smallest possible in order to
contain that data and avoid a large area spanning from -180 to 180
longitude. This means the resulting AreaDefinition will have a right-most
X extent greater than 180 degrees. This does not apply to data crossing
the north or south pole as there is no "smallest" area in this case.
Attributes:
area_id:
The name of the area.
description:
The description of the area.
projection:
The dictionary or string or CRS object of projection parameters.
Doesn't have to be complete. If not complete, ``proj_info`` must
be provided to ``freeze`` to "fill in" any missing parameters.
width:
x dimension in number of pixels, aka number of grid columns
height:
y dimension in number of pixels, aka number of grid rows
shape:
Corresponding array shape as (height, width)
area_extent:
The area extent of the area.
pixel_size_x:
Pixel width in projection units
pixel_size_y:
Pixel height in projection units
resolution:
Resolution of the resulting area as (pixel_size_x, pixel_size_y)
or a scalar if pixel_size_x == pixel_size_y.
optimize_projection:
Whether the projection parameters have to be optimized.
rotation:
Rotation in degrees (negative is cw)
"""

def __init__(self, area_id=None, description=None, projection=None,
Expand Down Expand Up @@ -958,21 +967,39 @@ def freeze(self, lonslats=None, resolution=None, shape=None, proj_info=None):
shape = None if None in shape else shape
area_extent = self.area_extent
if not area_extent or not width or not height:
proj4 = Proj(proj_dict)
try:
lons, lats = lonslats
except (TypeError, ValueError):
lons, lats = lonslats.get_lonlats()
xarr, yarr = proj4(np.asarray(lons), np.asarray(lats))
xarr[xarr > 9e29] = np.nan
yarr[yarr > 9e29] = np.nan
corners = [np.nanmin(xarr), np.nanmin(yarr),
np.nanmax(xarr), np.nanmax(yarr)]
corners = self._compute_bound_centers(proj_dict, lonslats)
area_extent, width, height = self.compute_domain(corners, resolution, shape)
return AreaDefinition(self.area_id, self.description, '',
projection, width, height,
area_extent, self.rotation)

def _compute_bound_centers(self, proj_dict, lonslats):
proj4 = Proj(proj_dict)
lons, lats = self._extract_lons_lats(lonslats)
# TODO: Do more dask-friendly things here
xarr, yarr = proj4(np.asarray(lons), np.asarray(lats))
xarr[xarr > 9e29] = np.nan
yarr[yarr > 9e29] = np.nan
xmin = np.nanmin(xarr)
xmax = np.nanmax(xarr)
ymin = np.nanmin(yarr)
ymax = np.nanmax(yarr)
x_passes_antimeridian = (xmax - xmin) > 355
epsilon = 0.1
y_is_pole = (ymax >= 90 - epsilon) or (ymin <= -90 + epsilon)
if proj4.crs.is_geographic and x_passes_antimeridian and not y_is_pole:
# cross anti-meridian of projection
xmin = np.nanmin(xarr[xarr >= 0])
xmax = np.nanmax(xarr[xarr < 0]) + 360
return xmin, ymin, xmax, ymax

def _extract_lons_lats(self, lonslats):
try:
lons, lats = lonslats
except (TypeError, ValueError):
lons, lats = lonslats.get_lonlats()
return lons, lats


def invproj(data_x, data_y, proj_dict):
"""Perform inverse projection."""
Expand Down
72 changes: 51 additions & 21 deletions pyresample/test/test_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -2291,7 +2291,7 @@ def _check_final_area_lon_lat_with_chunks(final_area, lons, lats, chunks):
np.testing.assert_array_equal(lats, lats_c)


class TestDynamicAreaDefinition(unittest.TestCase):
class TestDynamicAreaDefinition:
"""Test the DynamicAreaDefinition class."""

def test_freeze(self):
Expand All @@ -2308,10 +2308,10 @@ def test_freeze(self):
-872594.690447,
432079.38952,
904633.303964))
self.assertEqual(result.proj_dict['lon_0'], 16)
self.assertEqual(result.proj_dict['lat_0'], 58)
self.assertEqual(result.width, 288)
self.assertEqual(result.height, 592)
assert result.proj_dict['lon_0'] == 16
assert result.proj_dict['lat_0'] == 58
assert result.width == 288
assert result.height == 592

# make sure that setting `proj_info` once doesn't
# set it in the dynamic area
Expand All @@ -2322,11 +2322,41 @@ def test_freeze(self):
5380808.879250369,
1724415.6519203288,
6998895.701001488))
self.assertEqual(result.proj_dict['lon_0'], 0)
assert result.proj_dict['lon_0'] == 0
# lat_0 could be provided or not depending on version of pyproj
self.assertEqual(result.proj_dict.get('lat_0', 0), 0)
self.assertEqual(result.width, 395)
self.assertEqual(result.height, 539)
assert result.proj_dict.get('lat_0', 0) == 0
assert result.width == 395
assert result.height == 539

@pytest.mark.parametrize(
('lats',),
[
(np.linspace(-25.0, -10.0, 10),),
(np.linspace(10.0, 25.0, 10),),
(np.linspace(75, 90.0, 10),),
(np.linspace(-75, -90.0, 10),),
],
)
def test_freeze_longlat_antimeridian(self, lats):
"""Test geographic areas over the antimeridian."""
area = geometry.DynamicAreaDefinition('test_area', 'A test area',
'EPSG:4326')
lons = np.linspace(175, 185, 10)
lons[lons > 180] -= 360
result = area.freeze((lons, lats),
resolution=0.0056)

is_pole = (np.abs(lats) > 88).any()
extent = result.area_extent
if is_pole:
assert extent[0] < -178
assert extent[2] > 178
assert result.width == 64088
else:
assert extent[0] > 0
assert extent[2] > 0
assert result.width == 1787
assert result.height == 2680

def test_freeze_with_bb(self):
"""Test freezing the area with bounding box computation."""
Expand All @@ -2345,36 +2375,36 @@ def test_freeze_with_bb(self):
[-335439.956533, 5502125.451125,
191991.313351, 7737532.343683])

self.assertEqual(result.width, 4)
self.assertEqual(result.height, 18)
assert result.width == 4
assert result.height == 18
# Test for properties and shape usage in freeze.
area = geometry.DynamicAreaDefinition('test_area', 'A test area', {'proj': 'merc'},
width=4, height=18)
self.assertEqual((18, 4), area.shape)
assert (18, 4) == area.shape
result = area.freeze(sdef)
np.testing.assert_allclose(result.area_extent,
(996309.4426, 6287132.757981, 1931393.165263, 10837238.860543))
area = geometry.DynamicAreaDefinition('test_area', 'A test area', {'proj': 'merc'},
resolution=1000)
self.assertEqual(1000, area.pixel_size_x)
self.assertEqual(1000, area.pixel_size_y)
assert 1000 == area.pixel_size_x
assert 1000 == area.pixel_size_y

def test_compute_domain(self):
"""Test computing size and area extent."""
area = geometry.DynamicAreaDefinition('test_area', 'A test area',
{'proj': 'laea'})
corners = [1, 1, 9, 9]
self.assertRaises(ValueError, area.compute_domain, corners, 1, 1)
pytest.raises(ValueError, area.compute_domain, corners, 1, 1)

area_extent, x_size, y_size = area.compute_domain(corners, shape=(5, 5))
self.assertTupleEqual(area_extent, (0, 0, 10, 10))
self.assertEqual(x_size, 5)
self.assertEqual(y_size, 5)
assert area_extent == (0, 0, 10, 10)
assert x_size == 5
assert y_size == 5

area_extent, x_size, y_size = area.compute_domain(corners, resolution=2)
self.assertTupleEqual(area_extent, (0, 0, 10, 10))
self.assertEqual(x_size, 5)
self.assertEqual(y_size, 5)
assert area_extent == (0, 0, 10, 10)
assert x_size == 5
assert y_size == 5


class TestCrop(unittest.TestCase):
Expand Down

0 comments on commit 9801978

Please sign in to comment.