[Git][debian-gis-team/pyresample][master] 5 commits: New upstream version 1.23.0
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Sat Mar 26 06:54:53 GMT 2022
Antonio Valentino pushed to branch master at Debian GIS Project / pyresample
Commits:
5d1c4a63 by Antonio Valentino at 2022-03-26T06:12:40+00:00
New upstream version 1.23.0
- - - - -
7d55bfb8 by Antonio Valentino at 2022-03-26T06:13:16+00:00
Update upstream source from tag 'upstream/1.23.0'
Update to upstream version '1.23.0'
with Debian dir 92e17a3cc01dffa3d5a2d7f0cb8cebfa5be5c399
- - - - -
5b968484 by Antonio Valentino at 2022-03-26T06:14:27+00:00
New upstream release
- - - - -
dbef8cf5 by Antonio Valentino at 2022-03-26T06:17:42+00:00
Update d/copyright
- - - - -
2afad439 by Antonio Valentino at 2022-03-26T06:18:40+00:00
Set distribution to unstable
- - - - -
24 changed files:
- .github/workflows/ci.yaml
- CHANGELOG.md
- RELEASING.md
- debian/changelog
- debian/copyright
- pyresample/area_config.py
- pyresample/bilinear/_base.py
- pyresample/bilinear/xarr.py
- pyresample/ewa/_fornav_templates.cpp
- pyresample/geometry.py
- pyresample/gradient/__init__.py
- pyresample/spherical.py
- pyresample/spherical_geometry.py
- pyresample/test/test_bilinear.py
- pyresample/test/test_dask_ewa.py
- pyresample/test/test_files/areas.cfg
- pyresample/test/test_files/areas.yaml
- pyresample/test/test_geometry.py
- pyresample/test/test_gradient.py
- pyresample/test/test_spherical.py
- pyresample/test/test_spherical_geometry.py
- pyresample/test/test_utils.py
- pyresample/utils/rasterio.py
- pyresample/version.py
Changes:
=====================================
.github/workflows/ci.yaml
=====================================
@@ -46,7 +46,7 @@ jobs:
fail-fast: true
matrix:
os: ["windows-latest", "ubuntu-latest", "macos-latest"]
- python-version: ["3.7", "3.8", "3.9"]
+ python-version: ["3.8", "3.9", "3.10"]
experimental: [false]
include:
- python-version: "3.9"
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,65 @@
+## Version 1.23.0 (2022/03/21)
+
+### Issues Closed
+
+* [Issue 425](https://github.com/pytroll/pyresample/issues/425) - Pyresample/geometry.py resampling error related to dask.
+* [Issue 422](https://github.com/pytroll/pyresample/issues/422) - Cannot resample with `bilinear` from lat/lon grid onto MSG full disk ([PR 423](https://github.com/pytroll/pyresample/pull/423) by [@pnuu](https://github.com/pnuu))
+* [Issue 416](https://github.com/pytroll/pyresample/issues/416) - Unexpected results resampling Lambert Conformal to PlateCarree: pyresample or cartopy problem?
+
+In this release 3 issues were closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 426](https://github.com/pytroll/pyresample/pull/426) - Fix EWA resampling not ignoring fill values with maximum_weight_mode
+* [PR 424](https://github.com/pytroll/pyresample/pull/424) - Fix DynamicAreaDefinition resolution handling for incomplete projection definitions
+* [PR 423](https://github.com/pytroll/pyresample/pull/423) - Fix bilinear resampling to areas with invalid coordinates ([422](https://github.com/pytroll/pyresample/issues/422))
+* [PR 421](https://github.com/pytroll/pyresample/pull/421) - Fix inplace modification occuring in Arc.intersections
+* [PR 414](https://github.com/pytroll/pyresample/pull/414) - Fix gradient search for single band data
+
+#### Features added
+
+* [PR 415](https://github.com/pytroll/pyresample/pull/415) - Update AreaDefinition equality to use pyproj CRS
+* [PR 406](https://github.com/pytroll/pyresample/pull/406) - Change tested Python versions to 3.8, 3.9 and 3.10
+
+#### Backward incompatible changes
+
+* [PR 415](https://github.com/pytroll/pyresample/pull/415) - Update AreaDefinition equality to use pyproj CRS
+
+In this release 8 pull requests were closed.
+
+
+## Version 1.22.4 (2022/03/21)
+
+### Issues Closed
+
+* [Issue 425](https://github.com/pytroll/pyresample/issues/425) - Pyresample/geometry.py resampling error related to dask.
+* [Issue 416](https://github.com/pytroll/pyresample/issues/416) - Unexpected results resampling Lambert Conformal to PlateCarree: pyresample or cartopy problem?
+
+In this release 2 issues were closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 426](https://github.com/pytroll/pyresample/pull/426) - Fix EWA resampling not ignoring fill values with maximum_weight_mode
+* [PR 424](https://github.com/pytroll/pyresample/pull/424) - Fix DynamicAreaDefinition resolution handling for incomplete projection definitions
+* [PR 421](https://github.com/pytroll/pyresample/pull/421) - Fix inplace modification occuring in Arc.intersections
+* [PR 414](https://github.com/pytroll/pyresample/pull/414) - Fix gradient search for single band data
+
+#### Features added
+
+* [PR 415](https://github.com/pytroll/pyresample/pull/415) - Update AreaDefinition equality to use pyproj CRS
+* [PR 406](https://github.com/pytroll/pyresample/pull/406) - Change tested Python versions to 3.8, 3.9 and 3.10
+
+#### Backward incompatible changes
+
+* [PR 415](https://github.com/pytroll/pyresample/pull/415) - Update AreaDefinition equality to use pyproj CRS
+
+In this release 7 pull requests were closed.
+
+
## Version 1.22.3 (2021/12/07)
### Issues Closed
=====================================
RELEASING.md
=====================================
@@ -6,7 +6,7 @@
4. run `loghub` and update the `CHANGELOG.md` file:
```
-loghub pytroll/pyresample --token $LOGHUB_GITHUB_TOKEN -st v0.8.0 -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes" -plg backwards-incompatibility "Backward incompatible changes" -plg refactor "Refactoring"
+loghub pytroll/pyresample --token $LOGHUB_GITHUB_TOKEN -st $(git tag --sort=-version:refname --list 'v*' | head -n 1) -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes" -plg backwards-incompatibility "Backward incompatible changes" -plg refactor "Refactoring"
```
This uses a `LOGHUB_GITHUB_TOKEN` environment variable. This must be created
=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+pyresample (1.23.0-1) unstable; urgency=medium
+
+ * New upstream release.
+ * Update d/copyright.
+
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Sat, 26 Mar 2022 06:18:07 +0000
+
pyresample (1.22.3-2) unstable; urgency=medium
* Fix egg-info management in d/rules.
=====================================
debian/copyright
=====================================
@@ -36,7 +36,7 @@ Copyright: 2016-2019 David Hoese <david.hoese at ssec.wisc.edu>
License: GPL-3+
Files: debian/*
-Copyright: 2013-2021, Antonio Valentino <antonio.valentino at tiscali.it>
+Copyright: 2013-2022, Antonio Valentino <antonio.valentino at tiscali.it>
License: GPL-3+
License: GPL-3+
=====================================
pyresample/area_config.py
=====================================
@@ -477,7 +477,8 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No
p = Proj(crs, preserve_units=True)
except (RuntimeError, CRSError):
# Assume that an invalid projection will be "fixed" by a dynamic area definition later
- return _make_area(area_id, description, proj_id, projection, shape, area_extent, **kwargs)
+ return _make_area(area_id, description, proj_id, projection, shape, area_extent,
+ resolution=resolution, **kwargs)
# If no units are provided, try to get units used in proj_dict. If still none are provided, use meters.
if units is None:
=====================================
pyresample/bilinear/_base.py
=====================================
@@ -143,13 +143,15 @@ class BilinearBase(object):
self.mask_slices = self._index_array >= self._source_geo_def.size
def _get_valid_output_indices(self):
- return ((self._target_lons >= -180) & (self._target_lons <= 180) &
- (self._target_lats <= 90) & (self._target_lats >= -90))
+ self._valid_output_indices = np.ravel(
+ (self._target_lons >= -180) & (self._target_lons <= 180) &
+ (self._target_lats <= 90) & (self._target_lats >= -90))
def _get_index_array(self):
+ self._get_valid_output_indices()
index_array = _query_no_distance(
self._target_lons, self._target_lats,
- self._get_valid_output_indices(), self._resample_kdtree,
+ self._valid_output_indices, self._resample_kdtree,
self._neighbours, self._epsilon,
self._radius_of_influence)
self._index_array = self._reduce_index_array(index_array)
@@ -170,7 +172,10 @@ class BilinearBase(object):
corner_points, out_x, out_y)
def _get_output_xy(self):
- return _get_output_xy(self._target_geo_def)
+ out_x, out_y = _get_output_xy(self._target_geo_def)
+ out_x = out_x[self._valid_output_indices]
+ out_y = out_y[self._valid_output_indices]
+ return out_x, out_y
def _get_input_xy(self):
return _get_input_xy(self._source_geo_def,
@@ -637,9 +642,8 @@ def _resample(corner_points, fractional_distances):
def _query_no_distance(target_lons, target_lats,
valid_output_index, kdtree, neighbours, epsilon, radius):
"""Query the kdtree. No distances are returned."""
- voir = np.ravel(valid_output_index)
- target_lons_valid = np.ravel(target_lons)[voir]
- target_lats_valid = np.ravel(target_lats)[voir]
+ target_lons_valid = np.ravel(target_lons)[valid_output_index]
+ target_lats_valid = np.ravel(target_lats)[valid_output_index]
_, index_array = kdtree.query(
lonlat2xyz(target_lons_valid, target_lats_valid),
=====================================
pyresample/bilinear/xarr.py
=====================================
@@ -88,7 +88,10 @@ class XArrayBilinearResampler(BilinearBase):
self._valid_input_index, self._index_array)
def _get_output_xy(self):
- return _get_output_xy(self._target_geo_def)
+ out_x, out_y = _get_output_xy(self._target_geo_def)
+ out_x = out_x[self._valid_output_indices]
+ out_y = out_y[self._valid_output_indices]
+ return out_x, out_y
def _limit_output_values_to_input(self, data, res, fill_value):
epsilon = 1e-6
@@ -102,6 +105,19 @@ class XArrayBilinearResampler(BilinearBase):
return da.where(np.isnan(res), fill_value, res)
def _reshape_to_target_area(self, res, ndim):
+ if ndim == 3:
+ dim_multiplier = res.shape[0]
+ else:
+ dim_multiplier = 1
+ res = da.reshape(res, (1, res.size))
+ if res.size != dim_multiplier * self._target_geo_def.size:
+ out = []
+ for i in range(dim_multiplier):
+ tmp = da.full(self._target_geo_def.size, np.nan)
+ tmp[self._valid_output_indices] = res[i, :]
+ out.append(tmp)
+ res = da.stack(out)
+
shp = self._target_geo_def.shape
if ndim == 3:
res = da.reshape(res, (res.shape[0], shp[0], shp[1]))
=====================================
pyresample/ewa/_fornav_templates.cpp
=====================================
@@ -292,13 +292,9 @@ int compute_ewa(size_t chan_count, int maximum_weight_mode,
for (chan = 0; chan < chan_count; chan+=1) {
this_val = ((images[chan])[swath_offset]);
if (maximum_weight_mode) {
- if (weight > grid_weights[chan][grid_offset]) {
+ if (weight > grid_weights[chan][grid_offset] & !((this_val == img_fill) || (__isnan(this_val)))) {
((grid_weights[chan])[grid_offset]) = weight;
- if ((this_val == img_fill) || (__isnan(this_val))) {
- ((grid_accums[chan])[grid_offset]) = (accum_type)NPY_NANF;
- } else {
- ((grid_accums[chan])[grid_offset]) = (accum_type)this_val;
- }
+ ((grid_accums[chan])[grid_offset]) = (accum_type)this_val;
}
} else {
if ((this_val != img_fill) && !(__isnan(this_val))) {
@@ -405,13 +401,9 @@ int compute_ewa_single(int maximum_weight_mode,
this_val = (image[swath_offset]);
if (maximum_weight_mode) {
- if (weight > grid_weight[grid_offset]) {
+ if (weight > grid_weight[grid_offset] & !((this_val == img_fill) || (__isnan(this_val)))) {
grid_weight[grid_offset] = weight;
- if ((this_val == img_fill) || (__isnan(this_val))) {
- grid_accum[grid_offset] = (accum_type)NPY_NANF;
- } else {
- grid_accum[grid_offset] = (accum_type)this_val;
- }
+ grid_accum[grid_offset] = (accum_type)this_val;
}
} else {
if ((this_val != img_fill) && !(__isnan(this_val))) {
=====================================
pyresample/geometry.py
=====================================
@@ -827,7 +827,7 @@ class SwathDefinition(CoordinateDefinition):
new_proj.update(proj_dict)
return new_proj
- def _compute_uniform_shape(self):
+ def _compute_uniform_shape(self, resolution=None):
"""Compute the height and width of a domain to have uniform resolution across dimensions."""
g = Geod(ellps='WGS84')
@@ -861,14 +861,17 @@ class SwathDefinition(CoordinateDefinition):
az1, az2, width2 = g.inv(leftlons[-1], leftlats[-1], rightlons[-1], rightlats[-1])
az1, az2, height = g.inv(middlelons[0], middlelats[0], middlelons[-1], middlelats[-1])
width = min(width1, width2)
- vresolution = height * 1.0 / self.lons.shape[0]
- hresolution = width * 1.0 / self.lons.shape[1]
- resolution = min(vresolution, hresolution)
+ if resolution is None:
+ vresolution = height * 1.0 / self.lons.shape[0]
+ hresolution = width * 1.0 / self.lons.shape[1]
+ resolution = (vresolution, hresolution)
+ if isinstance(resolution, (tuple, list)):
+ resolution = min(*resolution)
width = int(width * 1.1 / resolution)
height = int(height * 1.1 / resolution)
return height, width
- def compute_optimal_bb_area(self, proj_dict=None):
+ def compute_optimal_bb_area(self, proj_dict=None, resolution=None):
"""Compute the "best" bounding box area for this swath with `proj_dict`.
By default, the projection is Oblique Mercator (`omerc` in proj.4), in
@@ -884,7 +887,7 @@ class SwathDefinition(CoordinateDefinition):
projection = proj_dict.setdefault('proj', 'omerc')
area_id = projection + '_otf'
description = 'On-the-fly ' + projection + ' area'
- height, width = self._compute_uniform_shape()
+ height, width = self._compute_uniform_shape(resolution)
proj_dict = self.compute_bb_proj_params(proj_dict)
area = DynamicAreaDefinition(area_id, description, proj_dict)
@@ -1024,7 +1027,7 @@ class DynamicAreaDefinition(object):
proj_info:
complementing parameters to the projection info.
- Resolution and shape parameters are ignored if the instance is created
+ Shape parameters are ignored if the instance is created
with the `optimize_projection` flag set to True.
"""
proj_dict = self._get_proj_dict()
@@ -1035,7 +1038,7 @@ class DynamicAreaDefinition(object):
projection = proj_dict
if self.optimize_projection:
- return lonslats.compute_optimal_bb_area(proj_dict)
+ return lonslats.compute_optimal_bb_area(proj_dict, resolution=resolution or self.resolution)
if resolution is None:
resolution = self.resolution
if shape is None:
@@ -1817,7 +1820,7 @@ class AreaDefinition(_ProjectionDefinition):
def __eq__(self, other):
"""Test for equality."""
try:
- return ((self.proj_str == other.proj_str) and
+ return ((self.crs == other.crs) and
(self.shape == other.shape) and
(np.allclose(self.area_extent, other.area_extent)))
except AttributeError:
=====================================
pyresample/gradient/__init__.py
=====================================
@@ -271,6 +271,9 @@ class GradientSearchResampler(BaseResampler):
if fill_value is not None:
res = da.where(np.isnan(res), fill_value, res)
+ if res.ndim > len(data_dims):
+ res = res.squeeze()
+
res = xr.DataArray(res, dims=data_dims, coords=coords)
return res
@@ -400,6 +403,6 @@ def _concatenate_chunks(chunks):
prev_y = y
res.append(da.concatenate(col, axis=1))
- res = da.concatenate(res, axis=2).squeeze()
+ res = da.concatenate(res, axis=2)
return res
=====================================
pyresample/spherical.py
=====================================
@@ -26,6 +26,13 @@ import numpy as np
logger = logging.getLogger(__name__)
+EPSILON = 0.0000001
+
+
+def _unwrap_radians(val, mod=np.pi):
+ """Put *val* between -*mod* and *mod*."""
+ return (val + mod) % (2 * mod) - mod
+
class SCoordinate(object):
"""Spherical coordinates.
@@ -35,7 +42,7 @@ class SCoordinate(object):
"""
def __init__(self, lon, lat):
- self.lon = lon
+ self.lon = float(_unwrap_radians(lon))
self.lat = lat
def cross2cart(self, point):
@@ -171,14 +178,6 @@ class CCoordinate(object):
np.arcsin(self.cart[2]))
-EPSILON = 0.0000001
-
-
-def modpi(val, mod=np.pi):
- """Put *val* between -*mod* and *mod*."""
- return (val + mod) % (2 * mod) - mod
-
-
class Arc(object):
"""An arc of the great circle between two points."""
@@ -252,17 +251,23 @@ class Arc(object):
From http://williams.best.vwh.net/intersect.htm
"""
+ end_lon = self.end.lon
+ other_end_lon = other_arc.end.lon
+
if self.end.lon - self.start.lon > np.pi:
- self.end.lon -= 2 * np.pi
+ end_lon -= 2 * np.pi
if other_arc.end.lon - other_arc.start.lon > np.pi:
- other_arc.end.lon -= 2 * np.pi
+ other_end_lon -= 2 * np.pi
if self.end.lon - self.start.lon < -np.pi:
- self.end.lon += 2 * np.pi
+ end_lon += 2 * np.pi
if other_arc.end.lon - other_arc.start.lon < -np.pi:
- other_arc.end.lon += 2 * np.pi
+ other_end_lon += 2 * np.pi
- ea_ = self.start.cross2cart(self.end).normalize()
- eb_ = other_arc.start.cross2cart(other_arc.end).normalize()
+ end_point = SCoordinate(end_lon, self.end.lat)
+ other_end_point = SCoordinate(other_end_lon, other_arc.end.lat)
+
+ ea_ = self.start.cross2cart(end_point).normalize()
+ eb_ = other_arc.start.cross2cart(other_end_point).normalize()
cross = ea_.cross(eb_)
lat = np.arctan2(cross.cart[2],
@@ -270,7 +275,7 @@ class Arc(object):
lon = np.arctan2(cross.cart[1], cross.cart[0])
return (SCoordinate(lon, lat),
- SCoordinate(modpi(lon + np.pi), -lat))
+ SCoordinate(_unwrap_radians(lon + np.pi), -lat))
def intersects(self, other_arc):
"""Check if the current arc and the *other_arc* intersect.
@@ -356,7 +361,7 @@ class SphPolygon:
radius (optional, number): Radius of spherical planet.
"""
self.vertices = vertices.astype(np.float64, copy=False)
- self.lon = self.vertices[:, 0]
+ self.lon = _unwrap_radians(self.vertices[:, 0])
self.lat = self.vertices[:, 1]
self.radius = radius
self.cvertices = np.array([np.cos(self.lat) * np.cos(self.lon),
=====================================
pyresample/spherical_geometry.py
=====================================
@@ -34,7 +34,11 @@ EPSILON = 0.0000001
class Coordinate(object):
- """Point on earth in terms of lat and lon."""
+ """Point on earth in terms of lat and lon.
+
+ It expects lon,lat in degrees
+ But self.lat and self.lon are returned in radians !
+ """
lat = None
lon = None
@@ -236,25 +240,30 @@ class Arc(object):
def intersections(self, other_arc):
"""Get the two intersections of the greats circles defined by the current arc and *other_arc*."""
+ end_lon = self.end.lon
+ other_end_lon = other_arc.end.lon
+
if self.end.lon - self.start.lon > math.pi:
- self.end.lon -= 2 * math.pi
+ end_lon -= 2 * math.pi
if other_arc.end.lon - other_arc.start.lon > math.pi:
- other_arc.end.lon -= 2 * math.pi
+ other_end_lon -= 2 * math.pi
if self.end.lon - self.start.lon < -math.pi:
- self.end.lon += 2 * math.pi
+ end_lon += 2 * math.pi
if other_arc.end.lon - other_arc.start.lon < -math.pi:
- other_arc.end.lon += 2 * math.pi
+ other_end_lon += 2 * math.pi
+
+ end_point = Coordinate(math.degrees(modpi(end_lon)), math.degrees(self.end.lat))
+ other_end_point = Coordinate(math.degrees(modpi(other_end_lon)), math.degrees(other_arc.end.lat))
- ea_ = self.start.cross2cart(self.end).normalize()
- eb_ = other_arc.start.cross2cart(other_arc.end).normalize()
+ ea_ = self.start.cross2cart(end_point).normalize()
+ eb_ = other_arc.start.cross2cart(other_end_point).normalize()
cross = ea_.cross(eb_)
lat = math.atan2(cross.z__, math.sqrt(cross.x__ ** 2 + cross.y__ ** 2))
lon = math.atan2(-cross.y__, cross.x__)
return (Coordinate(math.degrees(lon), math.degrees(lat)),
- Coordinate(math.degrees(modpi(lon + math.pi)),
- math.degrees(-lat)))
+ Coordinate(math.degrees(modpi(lon + math.pi)), math.degrees(-lat)))
def intersects(self, other_arc):
"""Determine if this arc and another arc intersect.
=====================================
pyresample/test/test_bilinear.py
=====================================
@@ -1032,7 +1032,8 @@ class TestXarrayBilinear(unittest.TestCase):
kdtree = mock.MagicMock()
kdtree.query.return_value = (1, 2)
lons, lats = self.target_def.get_lonlats()
- voi = (lons >= -180) & (lons <= 180) & (lats <= 90) & (lats >= -90)
+ voi = np.ravel(
+ (lons >= -180) & (lons <= 180) & (lats <= 90) & (lats >= -90))
res = _query_no_distance(lons, lats, voi, kdtree, self._neighbours,
0., self.radius)
# Only the second value from the query is returned
@@ -1159,3 +1160,47 @@ def test_check_fill_value():
# float fill value + float dtype -> no change
assert _check_fill_value(3.3, np.float32)
+
+
+def test_target_has_invalid_coordinates():
+ """Test bilinear resampling to area that has invalid coordinates.
+
+ The area used here is in geos projection that has space pixels in the corners.
+ """
+ import dask.array as da
+ import xarray as xr
+
+ from pyresample.bilinear import XArrayBilinearResampler
+
+ # NumpyBilinearResampler
+ from pyresample.geometry import AreaDefinition, GridDefinition
+
+ geos_def = AreaDefinition('geos',
+ 'GEO area with space in corners',
+ 'geos',
+ {'proj': 'geos',
+ 'lon_0': '0.0',
+ 'a': '6378169.0',
+ 'b': '6356583.8',
+ 'h': '35785831.0'},
+ 640, 640,
+ [-5432229.931711678,
+ -5429229.528545862,
+ 5429229.528545862,
+ 5432229.931711678])
+ lats = np.linspace(-89, 89, 179)
+ lons = np.linspace(-179, 179, 359)
+ lats = np.repeat(lats[:, None], 359, axis=1)
+ lons = np.repeat(lons[None, :], 179, axis=0)
+ grid_def = GridDefinition(lons=lons, lats=lats)
+
+ data_xr = xr.DataArray(da.random.uniform(0, 1, lons.shape), dims=["y", "x"])
+
+ resampler = XArrayBilinearResampler(grid_def,
+ geos_def,
+ 500e3,
+ reduce_data=False)
+ res = resampler.resample(data_xr)
+ res = res.compute()
+ assert not np.all(np.isnan(res))
+ assert np.any(np.isnan(res))
=====================================
pyresample/test/test_dask_ewa.py
=====================================
@@ -53,6 +53,9 @@ def _get_test_array(input_shape, input_dtype, chunk_size):
chunks=chunk_size, dtype=input_dtype)
else:
data = da.random.random(input_shape, chunks=chunk_size).astype(input_dtype)
+ fill_value = 127 if np.issubdtype(input_dtype, np.integer) else np.nan
+ if data.ndim in (2, 3):
+ data[..., int(data.shape[-2]) * 0.7, :] = fill_value
return data
@@ -204,6 +207,7 @@ class TestDaskEWAResampler:
mock.patch.object(source_swath, 'get_lonlats', wraps=source_swath.get_lonlats) as get_lonlats:
resampler = resampler_class(source_swath, target_area)
new_data = resampler.resample(swath_data, rows_per_scan=rows_per_scan,
+ weight_delta_max=40,
maximum_weight_mode=maximum_weight_mode)
_data_attrs_coords_checks(new_data, output_shape, input_dtype, target_area,
'test', 'test')
@@ -215,6 +219,7 @@ class TestDaskEWAResampler:
# resample a different dataset and make sure cache is used
swath_data2 = _create_second_test_data(swath_data)
new_data = resampler.resample(swath_data2, rows_per_scan=rows_per_scan,
+ weight_delta_max=40,
maximum_weight_mode=maximum_weight_mode)
_data_attrs_coords_checks(new_data, output_shape, input_dtype, target_area,
'test2', 'test2')
@@ -230,7 +235,11 @@ class TestDaskEWAResampler:
# check how many valid pixels we have
band_mult = 3 if 'bands' in result.dims else 1
fill_mask = _fill_mask(result.values)
- assert np.count_nonzero(~fill_mask) == 468 * band_mult
+ # without NaNs:
+ # exp_valid = 13939 if rows_per_scan == 10 else 14029
+ # with NaNs but no fix:
+ exp_valid = 13817 if rows_per_scan == 10 else 13913
+ assert np.count_nonzero(~fill_mask) == exp_valid * band_mult
@pytest.mark.parametrize(
('input_chunks', 'input_shape', 'input_dims'),
@@ -293,6 +302,7 @@ class TestDaskEWAResampler:
resampler = DaskEWAResampler(source_swath, target_area)
new_data = resampler.resample(swath_data, rows_per_scan=10,
+ weight_delta_max=40,
maximum_weight_mode=maximum_weight_mode)
assert new_data.shape == output_shape
assert new_data.dtype == np.float32
@@ -300,7 +310,7 @@ class TestDaskEWAResampler:
# check how many valid pixels we have
band_mult = 3 if len(output_shape) == 3 else 1
- assert np.count_nonzero(~np.isnan(new_data)) == 468 * band_mult
+ assert np.count_nonzero(~np.isnan(new_data)) == 13817 * band_mult
@pytest.mark.parametrize(
('input_shape', 'input_dims', 'maximum_weight_mode'),
=====================================
pyresample/test/test_files/areas.cfg
=====================================
@@ -37,7 +37,7 @@ REGION: pc_world {
REGION: ortho {
NAME: Ortho globe
PCS_ID: ortho_globe
- PCS_DEF: proj=ortho, a=6370997.0, lon_0=40, lat_0=-40
+ PCS_DEF: proj=ortho, lon_0=40, lat_0=-40, ellps=sphere
XSIZE: 640
YSIZE: 480
AREA_EXTENT: (-10000000, -10000000, 10000000, 10000000)
=====================================
pyresample/test/test_files/areas.yaml
=====================================
@@ -219,6 +219,14 @@ test_latlong:
lower_left_xy: [-0.08115781021773638, 0.4038691889114878]
upper_right_xy: [0.08115781021773638, 0.5427973973702365]
+omerc_bb_1000:
+ description: Oblique Mercator Bounding Box for Polar Overpasses
+ projection:
+ ellps: sphere
+ proj: omerc
+ optimize_projection: True
+ resolution: 1000
+
test_dynamic_resolution:
description: Dynamic with resolution specified in meters
projection:
=====================================
pyresample/test/test_geometry.py
=====================================
@@ -24,7 +24,6 @@ from unittest.mock import MagicMock, patch
import dask.array as da
import numpy as np
-import pyproj
import pytest
import xarray as xr
from pyproj import CRS
@@ -1287,7 +1286,7 @@ class Test(unittest.TestCase):
proj_dict['proj'] = 'omerc'
proj_dict['lat_0'] = 50.0
proj_dict['alpha'] = proj_dict.pop('lat_ts')
- proj_dict['no_rot'] = ''
+ proj_dict['no_rot'] = True
area = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
proj_dict, 10, 10,
[-1370912.72, -909968.64, 1029087.28,
@@ -1872,6 +1871,25 @@ class TestSwathDefinition(unittest.TestCase):
assert_np_dict_allclose(res.proj_dict, proj_dict)
self.assertEqual(res.shape, (6, 3))
+ def test_compute_optimal_bb_with_resolution(self):
+ """Test computing the bb area while passing in the resolution."""
+ import xarray as xr
+ nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+ [80.84000396728516, 60.74200439453125, 34.08500289916992],
+ [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+ lats = xr.DataArray(nplats)
+ nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+ [79.11000061035156, 7.284000396728516, -5.107000350952148],
+ [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+ lons = xr.DataArray(nplons)
+
+ area = geometry.SwathDefinition(lons, lats)
+
+ res1000 = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'}, resolution=1000)
+ res10000 = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'}, resolution=10000)
+ assert res1000.shape[0] // 10 == res10000.shape[0]
+ assert res1000.shape[1] // 10 == res10000.shape[1]
+
def test_aggregation(self):
"""Test aggregation on SwathDefinitions."""
import dask.array as da
@@ -2128,18 +2146,17 @@ class TestStackedAreaDefinition:
area1.crs, area1.width, y_size, area_extent)
@staticmethod
- def _compare_area_defs(actual, expected):
- actual_dict = actual.proj_dict
- if 'EPSG' in actual_dict or 'init' in actual_dict:
- # Use formal definition of EPSG projections to make them comparable to the base definition
- proj_def = pyproj.Proj(actual_dict).definition_string().strip()
- actual = actual.copy(projection=proj_def)
-
- # Remove extra attributes from the formal definition
- # for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']:
- # actual.proj_dict.pop(key, None)
-
- assert actual == expected
+ def _compare_area_defs(actual, expected, use_proj4=False):
+ if use_proj4:
+ # some EPSG codes have a lot of extra metadata that makes the CRS
+ # unequal. Skip real area equality and use this as an approximation
+ actual_str = actual.crs.to_proj4()
+ expected_str = expected.crs.to_proj4()
+ assert actual_str == expected_str
+ assert actual.shape == expected.shape
+ np.allclose(actual.area_extent, expected.area_extent)
+ else:
+ assert actual == expected
@pytest.mark.parametrize(
'projection',
@@ -2206,7 +2223,7 @@ class TestStackedAreaDefinition:
pytest.raises(ValueError, cad, *args, **kwargs)
else:
area_def = cad(*args, **kwargs)
- self._compare_area_defs(area_def, base_def)
+ self._compare_area_defs(area_def, base_def, use_proj4="EPSG" in projection)
def test_create_area_def_extra_combinations(self):
"""Test extra combinations of create_area_def parameters."""
@@ -2345,6 +2362,49 @@ class TestDynamicAreaDefinition:
assert result.width == 395
assert result.height == 539
+ def test_freeze_when_area_is_optimized_and_has_a_resolution(self):
+ """Test freezing an optimized area with a resolution."""
+ import xarray as xr
+ nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+ [80.84000396728516, 60.74200439453125, 34.08500289916992],
+ [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+ lats = xr.DataArray(nplats)
+ nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+ [79.11000061035156, 7.284000396728516, -5.107000350952148],
+ [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+ lons = xr.DataArray(nplons)
+
+ swath = geometry.SwathDefinition(lons, lats)
+
+ area10km = geometry.DynamicAreaDefinition('test_area', 'A test area',
+ {'ellps': 'WGS84', 'proj': 'omerc'},
+ resolution=10000,
+ optimize_projection=True)
+
+ result10km = area10km.freeze(swath)
+ assert result10km.shape == (679, 330)
+
+ def test_freeze_when_area_is_optimized_and_a_resolution_is_provided(self):
+ """Test freezing an optimized area when provided a resolution."""
+ import xarray as xr
+ nplats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+ [80.84000396728516, 60.74200439453125, 34.08500289916992],
+ [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+ lats = xr.DataArray(nplats)
+ nplons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+ [79.11000061035156, 7.284000396728516, -5.107000350952148],
+ [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+ lons = xr.DataArray(nplons)
+
+ swath = geometry.SwathDefinition(lons, lats)
+
+ area10km = geometry.DynamicAreaDefinition('test_area', 'A test area',
+ {'ellps': 'WGS84', 'proj': 'omerc'},
+ optimize_projection=True)
+
+ result10km = area10km.freeze(swath, 10000)
+ assert result10km.shape == (679, 330)
+
@pytest.mark.parametrize(
('lats',),
[
@@ -2401,7 +2461,7 @@ class TestDynamicAreaDefinition:
[50, 51, 52, 53]]
import xarray as xr
sdef = geometry.SwathDefinition(xr.DataArray(lons), xr.DataArray(lats))
- result = area.freeze(sdef, resolution=1000)
+ result = area.freeze(sdef)
np.testing.assert_allclose(result.area_extent,
[-335439.956533, 5502125.451125,
191991.313351, 7737532.343683])
=====================================
pyresample/test/test_gradient.py
=====================================
@@ -228,6 +228,16 @@ class TestGradientResampler(unittest.TestCase):
assert np.allclose(res[1, :, :], 2.0)
assert np.allclose(res[2, :, :], 3.0)
+ def test_resample_area_to_area_3d_single_channel(self):
+ """Resample area to area, 3d with only a single band."""
+ data = xr.DataArray(da.ones((1, ) + self.src_area.shape,
+ dtype=np.float64),
+ dims=['bands', 'y', 'x'])
+ res = self.resampler.compute(
+ data, method='bil').compute(scheduler='single-threaded')
+ assert res.shape == (1, ) + self.dst_area.shape
+ assert np.allclose(res[0, :, :], 1.0)
+
def test_resample_swath_to_area_2d(self):
"""Resample swath to area, 2d."""
data = xr.DataArray(da.ones(self.src_swath.shape, dtype=np.float64),
@@ -462,11 +472,11 @@ def test_concatenate_chunks():
(1, 1): [np.full((1, 3, 2), 0.5)],
(0, 1): [np.full((1, 3, 4), -1)]}
res = _concatenate_chunks(chunks).compute(scheduler='single-threaded')
- assert np.all(res[:5, :4] == 1.0)
- assert np.all(res[:5, 4:] == 0.0)
- assert np.all(res[5:, :4] == -1.0)
- assert np.all(res[5:, 4:] == 0.5)
- assert res.shape == (8, 6)
+ assert np.all(res[0, :5, :4] == 1.0)
+ assert np.all(res[0, :5, 4:] == 0.0)
+ assert np.all(res[0, 5:, :4] == -1.0)
+ assert np.all(res[0, 5:, 4:] == 0.5)
+ assert res.shape == (1, 8, 6)
# 3-band image
chunks = {(0, 0): [np.ones((3, 5, 4)), np.zeros((3, 5, 4))],
@@ -493,5 +503,4 @@ def test_concatenate_chunks_stack_calls(dask_da):
_ = _concatenate_chunks(chunks)
dask_da.stack.assert_called_once_with(chunks[(0, 0)], axis=-1)
dask_da.nanmax.assert_called_once()
- assert 'axis=2' in str(dask_da.concatenate.mock_calls[-2])
- assert 'squeeze' in str(dask_da.concatenate.mock_calls[-1])
+ assert 'axis=2' in str(dask_da.concatenate.mock_calls[-1])
=====================================
pyresample/test/test_spherical.py
=====================================
@@ -106,18 +106,29 @@ VERTICES_TEST_IS_INSIDE2 = np.array([[49.94506701, 46.52610743],
[49.94506701, 45.78080831]], dtype='float64')
+class TestRadiansFunctions(unittest.TestCase):
+ """Test utils functions processing radians arrays."""
+
+ def test_radian_unwrapping(self):
+ from pyresample.spherical import _unwrap_radians
+
+ # Test points at np.pi are converted to -np.pi
+ assert np.allclose(_unwrap_radians(-np.pi), -np.pi)
+ assert np.allclose(_unwrap_radians(np.pi), -np.pi)
+
+
class TestSCoordinate(unittest.TestCase):
"""Test SCoordinates."""
def test_distance(self):
"""Test Vincenty formula."""
d = SCoordinate(0, 0).distance(SCoordinate(1, 1))
- np.testing.assert_equal(d, 1.2745557823062943)
+ np.testing.assert_almost_equal(d, 1.2745557823062943)
def test_hdistance(self):
"""Test Haversine formula."""
d = SCoordinate(0, 0).hdistance(SCoordinate(1, 1))
- self.assertTrue(np.allclose(d, 1.2745557823062943))
+ np.testing.assert_allclose(d, 1.2745557823062943)
def test_str(self):
"""Check the string representation."""
@@ -129,6 +140,12 @@ class TestSCoordinate(unittest.TestCase):
d = SCoordinate(0, 0)
self.assertEqual(repr(d), "(0.0, 0.0)")
+ def test_equality_at_antimeridian(self):
+ """Test equality of points at the antimeridian."""
+ point_start = SCoordinate(-np.pi, 0)
+ point_end = SCoordinate(np.pi, 0)
+ assert point_start == point_end
+
class TestCCoordinate(unittest.TestCase):
"""Test SCoordinates."""
@@ -151,13 +168,13 @@ class TestCCoordinate(unittest.TestCase):
def test_normalize(self):
"""Normalize a cartesian vector."""
d = CCoordinate((2., 0., 0.))
- self.assertTrue(np.allclose(d.normalize().cart, [1, 0, 0]))
+ np.testing.assert_allclose(d.normalize().cart, [1, 0, 0])
def test_cross(self):
"""Test cross product in cartesian coordinates."""
d = CCoordinate((1., 0., 0.))
c = CCoordinate((0., 1., 0.))
- self.assertTrue(np.allclose(d.cross(c).cart, [0., 0., 1.]))
+ np.testing.assert_allclose(d.cross(c).cart, [0., 0., 1.])
def test_dot(self):
"""Test the dot product of two cartesian vectors."""
@@ -169,34 +186,31 @@ class TestCCoordinate(unittest.TestCase):
"""Test inequality of two cartesian vectors."""
d = CCoordinate((1., 0., 0.))
c = CCoordinate((0., 1., 0.))
- self.assertTrue(c != d)
+ assert c != d
def test_eq(self):
"""Test equality of two cartesian vectors."""
d = CCoordinate((1., 0., 0.))
c = CCoordinate((0., 1., 0.))
- self.assertFalse(c == d)
+ assert not c.__eq__(d)
def test_add(self):
"""Test adding cartesian vectors."""
d = CCoordinate((1., 0., 0.))
c = CCoordinate((0., 1., 0.))
b = CCoordinate((1., 1., 0.))
- self.assertTrue(np.allclose((d + c).cart, b.cart))
-
- self.assertTrue(np.allclose((d + (0, 1, 0)).cart, b.cart))
-
- self.assertTrue(np.allclose(((0, 1, 0) + d).cart, b.cart))
+ np.testing.assert_allclose((d + c).cart, b.cart)
+ np.testing.assert_allclose((d + (0, 1, 0)).cart, b.cart)
+ np.testing.assert_allclose(((0, 1, 0) + d).cart, b.cart)
def test_mul(self):
"""Test multiplying (element-wise) cartesian vectors."""
d = CCoordinate((1., 0., 0.))
c = CCoordinate((0., 1., 0.))
b = CCoordinate((0., 0., 0.))
- self.assertTrue(np.allclose((d * c).cart, b.cart))
- self.assertTrue(np.allclose((d * (0, 1, 0)).cart, b.cart))
-
- self.assertTrue(np.allclose(((0, 1, 0) * d).cart, b.cart))
+ np.testing.assert_allclose((d * c).cart, b.cart)
+ np.testing.assert_allclose((d * (0, 1, 0)).cart, b.cart)
+ np.testing.assert_allclose(((0, 1, 0) * d).cart, b.cart)
def test_to_spherical(self):
"""Test converting to spherical coordinates."""
@@ -214,9 +228,8 @@ class TestArc(unittest.TestCase):
arc2 = Arc(SCoordinate(0, np.deg2rad(10)),
SCoordinate(np.deg2rad(10), 0))
- self.assertFalse(arc1 == arc2)
-
- self.assertTrue(arc1 == arc1)
+ assert not arc1.__eq__(arc2)
+ assert arc1 == arc1
def test_ne(self):
arc1 = Arc(SCoordinate(0, 0),
@@ -224,9 +237,8 @@ class TestArc(unittest.TestCase):
arc2 = Arc(SCoordinate(0, np.deg2rad(10)),
SCoordinate(np.deg2rad(10), 0))
- self.assertTrue(arc1 != arc2)
-
- self.assertFalse(arc1 != arc1)
+ assert arc1 != arc2
+ assert not arc1.__ne__(arc1)
def test_str(self):
arc1 = Arc(SCoordinate(0, 0),
@@ -241,13 +253,13 @@ class TestArc(unittest.TestCase):
SCoordinate(np.deg2rad(10), 0))
lon, lat = arc1.intersection(arc2)
- self.assertTrue(np.allclose(np.rad2deg(lon), 5))
+ np.testing.assert_allclose(np.rad2deg(lon), 5)
self.assertEqual(np.rad2deg(lat).round(7), round(5.0575148968282093, 7))
arc1 = Arc(SCoordinate(0, 0),
SCoordinate(np.deg2rad(10), np.deg2rad(10)))
- self.assertTrue(arc1.intersection(arc1) is None)
+ assert (arc1.intersection(arc1) is None)
arc1 = Arc(SCoordinate(np.deg2rad(24.341215776575297),
np.deg2rad(44.987819588259327)),
@@ -270,7 +282,7 @@ class TestArc(unittest.TestCase):
SCoordinate(np.deg2rad(-5.893976312685715),
np.deg2rad(48.445795283217116)))
- self.assertTrue(arc1.intersection(arc2) is None)
+ assert (arc1.intersection(arc2) is None)
def test_angle(self):
arc1 = Arc(SCoordinate(np.deg2rad(157.5),
@@ -308,6 +320,40 @@ class TestArc(unittest.TestCase):
arc2 = Arc(SCoordinate(2, 0), SCoordinate(3, 0))
self.assertRaises(ValueError, arc1.angle, arc2)
+ def test_disjoint_arcs_crossing_antimeridian(self):
+ import copy
+ arc1 = Arc(SCoordinate(*np.deg2rad((143.76, 0))),
+ SCoordinate(*np.deg2rad((143.95, 7.33)))
+ )
+ arc2 = Arc(SCoordinate(*np.deg2rad((170.34, 71.36))),
+ SCoordinate(*np.deg2rad((-171.03, 76.75)))
+ )
+ arc1_orig = copy.deepcopy(arc1)
+ arc2_orig = copy.deepcopy(arc2)
+ point = arc1.intersection(arc2)
+ # Assert original arcs are unaffected
+ assert arc1_orig.end.lon == arc1.end.lon
+ assert arc2_orig.end.lon == arc2.end.lon
+ # Assert disjoint arcs returns None
+ assert isinstance(point, type(None))
+
+ def test_intersecting_arcs_crossing_antimeridian(self):
+ import copy
+ arc1 = Arc(SCoordinate(*np.deg2rad((-180.0, -90.0))),
+ SCoordinate(*np.deg2rad((-180.0, 90.0)))
+ )
+ arc2 = Arc(SCoordinate(*np.deg2rad((-171.03, -76.75))),
+ SCoordinate(*np.deg2rad((170.34, -71.36)))
+ )
+ arc1_orig = copy.deepcopy(arc1)
+ arc2_orig = copy.deepcopy(arc2)
+ point = arc1.intersection(arc2)
+ # Assert original arcs are unaffected
+ assert arc1_orig.end.lon == arc1.end.lon
+ assert arc2_orig.end.lon == arc2.end.lon
+ # Assert intersection result
+ assert point == SCoordinate(*np.deg2rad((-180.0, -74.7884716)))
+
class TestSphericalPolygon(unittest.TestCase):
"""Test the spherical polygon."""
@@ -383,41 +429,41 @@ class TestSphericalPolygon(unittest.TestCase):
polygon2 = SphPolygon(np.deg2rad(vertices))
- self.assertTrue(polygon1._is_inside(polygon2))
+ assert polygon1._is_inside(polygon2)
- self.assertFalse(polygon2._is_inside(polygon1))
+ assert not polygon2._is_inside(polygon1)
# Why checking the areas here!? It has nothing to do with the is_inside function!
- self.assertTrue(polygon2.area() > polygon1.area())
+ assert polygon2.area() > polygon1.area()
polygon2.invert()
- self.assertFalse(polygon1._is_inside(polygon2))
- self.assertFalse(polygon2._is_inside(polygon1))
+ assert not polygon1._is_inside(polygon2)
+ assert not polygon2._is_inside(polygon1)
vertices = np.array([[0, 0, 30, 30],
[21, 30, 30, 21]]).T
polygon2 = SphPolygon(np.deg2rad(vertices))
- self.assertFalse(polygon1._is_inside(polygon2))
- self.assertFalse(polygon2._is_inside(polygon1))
+ assert not polygon1._is_inside(polygon2)
+ assert not polygon2._is_inside(polygon1)
polygon2.invert()
- self.assertTrue(polygon1._is_inside(polygon2))
- self.assertFalse(polygon2._is_inside(polygon1))
+ assert polygon1._is_inside(polygon2)
+ assert not polygon2._is_inside(polygon1)
vertices = np.array([[100, 100, 130, 130],
[41, 50, 50, 41]]).T
polygon2 = SphPolygon(np.deg2rad(vertices))
- self.assertFalse(polygon1._is_inside(polygon2))
- self.assertFalse(polygon2._is_inside(polygon1))
+ assert not polygon1._is_inside(polygon2)
+ assert not polygon2._is_inside(polygon1)
polygon2.invert()
- self.assertTrue(polygon1._is_inside(polygon2))
- self.assertFalse(polygon2._is_inside(polygon1))
+ assert polygon1._is_inside(polygon2)
+ assert not polygon2._is_inside(polygon1)
vertices = VERTICES_TEST_IS_INSIDE1
@@ -427,8 +473,8 @@ class TestSphericalPolygon(unittest.TestCase):
polygon2 = SphPolygon(np.deg2rad(vertices))
- self.assertFalse(polygon2._is_inside(polygon1))
- self.assertFalse(polygon1._is_inside(polygon2))
+ assert not polygon2._is_inside(polygon1)
+ assert not polygon1._is_inside(polygon2)
def test_is_inside_float32(self):
"""Test that precision dependent calculations work.
@@ -478,7 +524,7 @@ class TestSphericalPolygon(unittest.TestCase):
poly2 = SphPolygon(np.deg2rad(vertices))
uni = np.array([[157.5, 89.23460094],
- [-225., 89.],
+ [135., 89.],
[112.5, 89.23460094],
[90., 89.],
[67.5, 89.23460094],
@@ -495,8 +541,7 @@ class TestSphericalPolygon(unittest.TestCase):
[-180., 89.]])
poly_union = poly1.union(poly2)
-
- self.assertTrue(np.allclose(poly_union.vertices, np.deg2rad(uni)))
+ assert np.allclose(poly_union.vertices, np.deg2rad(uni))
def test_union_polygons_overlaps_completely(self):
"""Test the union method when one polygon is entirely inside the other."""
@@ -515,9 +560,8 @@ class TestSphericalPolygon(unittest.TestCase):
expected = np.deg2rad(np.array([[0, 0, 30, 30],
[0, 30, 30, 0]]).T)
- self.assertTrue(np.allclose(poly_union1.vertices, expected))
-
- self.assertTrue(np.allclose(poly_union2.vertices, expected))
+ np.testing.assert_allclose(poly_union1.vertices, expected)
+ np.testing.assert_allclose(poly_union2.vertices, expected)
def test_intersection(self):
"""Test the intersection function."""
@@ -538,8 +582,7 @@ class TestSphericalPolygon(unittest.TestCase):
[-157.5, 89.23460094]])
poly_inter = poly1.intersection(poly2)
- self.assertTrue(np.allclose(poly_inter.vertices,
- np.deg2rad(inter)))
+ np.testing.assert_allclose(poly_inter.vertices, np.deg2rad(inter))
# Test 2 polygons sharing 2 contiguous edges.
@@ -567,8 +610,7 @@ class TestSphericalPolygon(unittest.TestCase):
poly2 = SphPolygon(np.deg2rad(vertices2))
poly_inter = poly1.intersection(poly2)
- self.assertTrue(np.allclose(poly_inter.vertices,
- np.deg2rad(vertices3)))
+ np.testing.assert_allclose(poly_inter.vertices, np.deg2rad(vertices3))
# Test when last node of the intersection is the last vertice of the
# second polygon.
@@ -598,12 +640,38 @@ class TestSphericalPolygon(unittest.TestCase):
poly2 = SphPolygon(np.deg2rad(area_vertices))
poly_inter = poly1.intersection(poly2)
- self.assertTrue(np.allclose(poly_inter.vertices,
- np.deg2rad(res)))
+ np.testing.assert_allclose(poly_inter.vertices, np.deg2rad(res))
poly_inter = poly2.intersection(poly1)
- self.assertTrue(np.allclose(poly_inter.vertices,
- np.deg2rad(res)))
+ np.testing.assert_allclose(poly_inter.vertices, np.deg2rad(res))
+
+ # Test when intersection occurs across the antimeridian
+ vertices_epsg_4326 = np.deg2rad(np.array([[-180, -90],
+ [-180, 90],
+ [0, 90],
+ [180, 90],
+ [180, -90]]))
+ vertices_goes_west = np.array([[2.50919361e+00, 1.19782309e-16],
+ [2.51252770e+00, 1.27991660e-01],
+ [2.97300443e+00, 1.24550237e+00],
+ [-2.98515478e+00, 1.33966326e+00],
+ [-1.00821045e+00, -0.00000000e+00],
+ [-2.98515478e+00, -1.33966326e+00],
+ [2.97300443e+00, -1.24550237e+00],
+ [2.51252770e+00, -1.27991660e-01]])
+
+ goes_west_poly = SphPolygon(vertices_goes_west)
+ epsg_4326_poly = SphPolygon(vertices_epsg_4326)
+
+ intersection_poly = epsg_4326_poly.intersection(goes_west_poly)
+
+ # Assert found the correct intersection point
+ assert np.allclose(intersection_poly.vertices[2, :],
+ np.array([-2.98515478, 1.33966326]))
+
+ # Check bounded between -180 and 180
+ assert np.all(intersection_poly.vertices >= -np.pi)
+ assert np.all(intersection_poly.vertices <= np.pi)
def test_consistent_radius(self):
poly1 = np.array([(-50, 69), (-36, 69), (-36, 64), (-50, 64)])
=====================================
pyresample/test/test_spherical_geometry.py
=====================================
@@ -45,23 +45,23 @@ class TestOverlap(unittest.TestCase):
point = Coordinate(0, 0)
- self.assertTrue(point in area)
+ assert (point in area)
point = Coordinate(0, 12)
- self.assertFalse(point in area)
+ assert not (point in area)
lons = np.array([[-179, 179], [-179, 179]])
lats = np.array([[1, 1], [-1, -1]])
area = geometry.SwathDefinition(lons, lats)
point = Coordinate(180, 0)
- self.assertTrue(point in area)
+ assert (point in area)
point = Coordinate(180, 12)
- self.assertFalse(point in area)
+ assert not (point in area)
point = Coordinate(-180, 12)
- self.assertFalse(point in area)
+ assert not (point in area)
self.assert_raises(ValueError, Coordinate, 0, 192)
@@ -73,7 +73,7 @@ class TestOverlap(unittest.TestCase):
area = geometry.SwathDefinition(lons, lats)
point = Coordinate(90, 90)
- self.assertTrue(point in area)
+ assert (point in area)
def test_overlaps(self):
"""Test if two areas overlap."""
@@ -85,8 +85,8 @@ class TestOverlap(unittest.TestCase):
lats2 = np.array([[89, 89], [89, 89]])
area2 = geometry.SwathDefinition(lons2, lats2)
- self.assertTrue(area1.overlaps(area2))
- self.assertTrue(area2.overlaps(area1))
+ assert area1.overlaps(area2)
+ assert area2.overlaps(area1)
lons1 = np.array([[0, 45], [135, 90]])
lats1 = np.array([[89, 89], [89, 89]])
@@ -96,8 +96,8 @@ class TestOverlap(unittest.TestCase):
lats2 = np.array([[89, 89], [89, 89]])
area2 = geometry.SwathDefinition(lons2, lats2)
- self.assertFalse(area1.overlaps(area2))
- self.assertFalse(area2.overlaps(area1))
+ assert not area1.overlaps(area2)
+ assert not area2.overlaps(area1)
lons1 = np.array([[-1, 1], [-1, 1]])
lats1 = np.array([[1, 1], [-1, -1]])
@@ -107,8 +107,8 @@ class TestOverlap(unittest.TestCase):
lats2 = np.array([[0, 0], [2, 2]])
area2 = geometry.SwathDefinition(lons2, lats2)
- self.assertTrue(area1.overlaps(area2))
- self.assertTrue(area2.overlaps(area1))
+ assert area1.overlaps(area2)
+ assert area2.overlaps(area1)
lons1 = np.array([[-1, 0], [-1, 0]])
lats1 = np.array([[1, 2], [-1, 0]])
@@ -118,8 +118,8 @@ class TestOverlap(unittest.TestCase):
lats2 = np.array([[1, 2], [-1, 0]])
area2 = geometry.SwathDefinition(lons2, lats2)
- self.assertFalse(area1.overlaps(area2))
- self.assertFalse(area2.overlaps(area1))
+ assert not area1.overlaps(area2)
+ assert not area2.overlaps(area1)
def test_overlap_rate(self):
"""Test how much two areas overlap."""
@@ -333,17 +333,17 @@ class TestSphereGeometry(unittest.TestCase):
arc35 = Arc(p3_, p5_)
- self.assertTrue(arc13.intersects(arc24))
+ assert arc13.intersects(arc24)
- self.assertFalse(arc32.intersects(arc41))
+ assert not arc32.intersects(arc41)
- self.assertFalse(arc56.intersects(arc40))
+ assert not arc56.intersects(arc40)
- self.assertFalse(arc56.intersects(arc40))
+ assert not arc56.intersects(arc40)
- self.assertFalse(arc45.intersects(arc02))
+ assert not arc45.intersects(arc02)
- self.assertTrue(arc35.intersects(arc24))
+ assert arc35.intersects(arc24)
p0_ = Coordinate(180, 0)
p1_ = Coordinate(180, 1)
@@ -367,17 +367,17 @@ class TestSphereGeometry(unittest.TestCase):
arc35 = Arc(p3_, p5_)
- self.assertTrue(arc13.intersects(arc24))
+ assert arc13.intersects(arc24)
- self.assertFalse(arc32.intersects(arc41))
+ assert not arc32.intersects(arc41)
- self.assertFalse(arc56.intersects(arc40))
+ assert not arc56.intersects(arc40)
- self.assertFalse(arc56.intersects(arc40))
+ assert not arc56.intersects(arc40)
- self.assertFalse(arc45.intersects(arc02))
+ assert not arc45.intersects(arc02)
- self.assertTrue(arc35.intersects(arc24))
+ assert arc35.intersects(arc24)
# case of the north pole
@@ -403,14 +403,14 @@ class TestSphereGeometry(unittest.TestCase):
arc35 = Arc(p3_, p5_)
- self.assertTrue(arc13.intersects(arc24))
+ assert arc13.intersects(arc24)
- self.assertFalse(arc32.intersects(arc41))
+ assert not arc32.intersects(arc41)
- self.assertFalse(arc56.intersects(arc40))
+ assert not arc56.intersects(arc40)
- self.assertFalse(arc56.intersects(arc40))
+ assert not arc56.intersects(arc40)
- self.assertFalse(arc45.intersects(arc02))
+ assert not arc45.intersects(arc02)
- self.assertTrue(arc35.intersects(arc24))
+ assert arc35.intersects(arc24)
=====================================
pyresample/test/test_utils.py
=====================================
@@ -191,6 +191,13 @@ Area extent: (-0.0812, 0.4039, 0.0812, 0.5428)""".format(projection)
self.assertTrue(hasattr(test_area, 'resolution'))
np.testing.assert_allclose(test_area.resolution, (1.0, 1.0))
+ def test_dynamic_area_parser_passes_resolution(self):
+ """Test that the resolution from the file is passed to a dynamic area."""
+ from pyresample import parse_area_file
+ test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area = parse_area_file(test_area_file, 'omerc_bb_1000')[0]
+ assert test_area.resolution == (1000, 1000)
+
def test_multiple_file_content(self):
from pyresample import parse_area_file
from pyresample.area_config import load_area_from_string
@@ -492,8 +499,7 @@ class TestMisc(unittest.TestCase):
transform = Affine(300.0379266750948, 0.0, 101985.0,
0.0, -300.041782729805, 2826915.0)
source = tmptiff(transform=transform)
- proj_dict = {'init': 'epsg:3857'}
- area_def = utils.rasterio.get_area_def_from_raster(source, proj_dict=proj_dict)
+ area_def = utils.rasterio.get_area_def_from_raster(source, projection="EPSG:3857")
self.assertEqual(area_def.crs, CRS(3857))
=====================================
pyresample/utils/rasterio.py
=====================================
@@ -19,7 +19,7 @@
from . import proj4_str_to_dict
-def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, proj_dict=None):
+def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, projection=None):
from pyresample.geometry import AreaDefinition
# a: width of a pixel
@@ -33,46 +33,44 @@ def _get_area_def_from_gdal(dataset, area_id=None, name=None, proj_id=None, proj
raise ValueError('Rotated rasters are not supported at this time.')
area_extent = (c, f + e * dataset.RasterYSize, c + a * dataset.RasterXSize, f)
- if proj_dict is None:
+ if projection is None:
from osgeo import osr
proj = dataset.GetProjection()
if proj != '':
sref = osr.SpatialReference(wkt=proj)
- proj_dict = proj4_str_to_dict(sref.ExportToProj4())
+ projection = proj4_str_to_dict(sref.ExportToProj4())
else:
raise ValueError('The source raster is not gereferenced, please provide the value of proj_dict')
if proj_id is None:
proj_id = proj.split('"')[1]
- area_def = AreaDefinition(area_id, name, proj_id, proj_dict,
+ area_def = AreaDefinition(area_id, name, proj_id, projection,
dataset.RasterXSize, dataset.RasterYSize, area_extent)
return area_def
-def _get_area_def_from_rasterio(dataset, area_id, name, proj_id=None, proj_dict=None):
+def _get_area_def_from_rasterio(dataset, area_id, name, proj_id=None, projection=None):
from pyresample.geometry import AreaDefinition
a, b, c, d, e, f, _, _, _ = dataset.transform
if not (b == d == 0):
raise ValueError('Rotated rasters are not supported at this time.')
- if proj_dict is None:
- crs = dataset.crs
- if crs is not None:
- proj_dict = dataset.crs.to_dict()
- else:
+ if projection is None:
+ projection = dataset.crs
+ if projection is None:
raise ValueError('The source raster is not gereferenced, please provide the value of proj_dict')
if proj_id is None:
- proj_id = crs.wkt.split('"')[1]
+ proj_id = projection.wkt.split('"')[1]
- area_def = AreaDefinition(area_id, name, proj_id, proj_dict,
+ area_def = AreaDefinition(area_id, name, proj_id, projection,
dataset.width, dataset.height, dataset.bounds)
return area_def
-def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj_dict=None):
+def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, projection=None):
"""Construct AreaDefinition object from raster.
Parameters
@@ -86,8 +84,9 @@ def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj
Name of area
proj_id : str, optional
ID of projection
- proj_dict : dict, optional
- PROJ.4 parameters
+ projection : pyproj.CRS, dict, or string, optional
+ Projection parameters as a CRS object or a dictionary or string of
+ PROJ.4 parameters.
Returns
-------
@@ -114,8 +113,8 @@ def get_area_def_from_raster(source, area_id=None, name=None, proj_id=None, proj
try:
if rasterio is not None and isinstance(source, (rasterio.io.DatasetReader, rasterio.io.DatasetWriter)):
- return _get_area_def_from_rasterio(source, area_id, name, proj_id, proj_dict)
- return _get_area_def_from_gdal(source, area_id, name, proj_id, proj_dict)
+ return _get_area_def_from_rasterio(source, area_id, name, proj_id, projection)
+ return _get_area_def_from_gdal(source, area_id, name, proj_id, projection)
finally:
if cleanup_rasterio:
source.close()
=====================================
pyresample/version.py
=====================================
@@ -23,9 +23,9 @@ def get_keywords():
# setup.py/versioneer.py will grep for the variable names, so they must
# each be defined on a line of their own. _version.py will just call
# get_keywords().
- git_refnames = " (HEAD -> main, tag: v1.22.3)"
- git_full = "3e8d2442790455493e6c376408e90607bd485aad"
- git_date = "2021-12-07 20:02:34 -0600"
+ git_refnames = " (tag: v1.23.0)"
+ git_full = "3a3d3a91a5e30445a9c5ce45d1c5ccedf07e9140"
+ git_date = "2022-03-21 20:31:11 -0500"
keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
return keywords
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyresample/-/compare/c5d6e139f687ec590dd3e81a6277f59aa0d42907...2afad43957077f743b4c91ef8b86afc54521b850
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyresample/-/compare/c5d6e139f687ec590dd3e81a6277f59aa0d42907...2afad43957077f743b4c91ef8b86afc54521b850
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20220326/13f84f52/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list