[Git][debian-gis-team/pyresample][upstream] New upstream version 1.29.0
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Sun Aug 4 17:15:55 BST 2024
Antonio Valentino pushed to branch upstream at Debian GIS Project / pyresample
Commits:
95525b45 by Antonio Valentino at 2024-08-04T15:24:03+00:00
New upstream version 1.29.0
- - - - -
11 changed files:
- .github/workflows/deploy.yaml
- .pre-commit-config.yaml
- CHANGELOG.md
- pyresample/_formatting_html.py
- pyresample/bucket/__init__.py
- pyresample/geometry.py
- pyresample/static/css/style.css
- pyresample/test/test_bucket.py
- pyresample/test/test_geometry/test_area.py
- pyresample/utils/cartopy.py
- pyresample/version.py
Changes:
=====================================
.github/workflows/deploy.yaml
=====================================
@@ -66,7 +66,7 @@ jobs:
platforms: all
- name: Build wheels
- uses: pypa/cibuildwheel at v2.19.1
+ uses: pypa/cibuildwheel at v2.19.2
env:
CIBW_SKIP: "cp36-* cp37-* cp38-* pp* *i686 *-musllinux*"
CIBW_ARCHS: "${{ matrix.cibw_archs }}"
=====================================
.pre-commit-config.yaml
=====================================
@@ -2,12 +2,12 @@ exclude: '^$'
fail_fast: false
repos:
- repo: https://github.com/astral-sh/ruff-pre-commit
- rev: 'v0.4.7'
+ rev: 'v0.5.0'
hooks:
- id: ruff
# once https://github.com/astral-sh/ruff/issues/2402 is fully resolved then we can get rid of flake8:
- repo: https://github.com/pycqa/flake8
- rev: 7.0.0
+ rev: 7.1.0
hooks:
- id: flake8
- repo: https://github.com/pre-commit/pre-commit-hooks
@@ -18,12 +18,12 @@ repos:
- id: check-yaml
args: [--unsafe]
- repo: https://github.com/PyCQA/bandit
- rev: '1.7.8' # Update me!
+ rev: '1.7.9' # Update me!
hooks:
- id: bandit
args: [--ini, .bandit]
- repo: https://github.com/pre-commit/mirrors-mypy
- rev: 'v1.10.0' # Use the sha / tag you want to point at
+ rev: 'v1.10.1' # Use the sha / tag you want to point at
hooks:
- id: mypy
additional_dependencies:
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,27 @@
+## Version 1.29.0 (2024/07/31)
+
+### Issues Closed
+
+* [Issue 609](https://github.com/pytroll/pyresample/issues/609) - Error in SwathDefinition html representation if lon/lat arrays are dask arrays ([PR 610](https://github.com/pytroll/pyresample/pull/610) by [@BENR0](https://github.com/BENR0))
+* [Issue 354](https://github.com/pytroll/pyresample/issues/354) - `get_sum` not matched with `bucket_sum`
+
+In this release 2 issues were closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 610](https://github.com/pytroll/pyresample/pull/610) - Fix SwathDefinition html representation error when lons/lats 1D ([609](https://github.com/pytroll/pyresample/issues/609))
+* [PR 601](https://github.com/pytroll/pyresample/pull/601) - Fix spill over of ocean/and features in cartopy plots in case of geostationary full disc plot.
+* [PR 596](https://github.com/pytroll/pyresample/pull/596) - Fix AreaDefinition array index methods mishandling outer edge values ([691](https://github.com/ssec/polar2grid/issues/691))
+
+#### Features added
+
+* [PR 602](https://github.com/pytroll/pyresample/pull/602) - Add support for `fill_value` and `set_empty_bucket_to` in BucketResampler `get_sum`
+
+In this release 4 pull requests were closed.
+
+
## Version 1.28.4 (2024/07/01)
### Pull Requests Merged
=====================================
pyresample/_formatting_html.py
=====================================
@@ -25,6 +25,7 @@ from typing import Literal, Optional, Union
import numpy as np
import pyresample.geometry as geom
+from pyresample.utils.proj4 import ignore_pyproj_proj_warnings
try:
import cartopy
@@ -65,15 +66,11 @@ def _icon(icon_name):
def plot_area_def(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'], # noqa F821
- feature_res: Optional[str] = "110m",
fmt: Optional[Literal["svg", "png", None]] = None) -> Union[str, None]:
"""Plot area.
Args:
area : Area/Swath to plot.
- feature_res :
- Resolution of the features added to the map. Argument is handed over
- to `scale` parameter in cartopy.feature.
fmt : Output format of the plot. The output is the string representation of
the respective format xml for svg and base64 for png. Either svg or png.
If None (default) plot is just shown.
@@ -102,21 +99,10 @@ def plot_area_def(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'], #
bounds = poly.buffer(5).bounds
ax.set_extent([bounds[0], bounds[2], bounds[1], bounds[3]], crs=cartopy.crs.CRS(area.crs))
- coastlines = cartopy.feature.NaturalEarthFeature(category="physical",
- name="coastline",
- scale=feature_res,
- linewidth=1,
- facecolor="never")
- borders = cartopy.feature.NaturalEarthFeature(category="cultural",
- name="admin_0_boundary_lines_land", # noqa E114
- scale=feature_res,
- edgecolor="black",
- facecolor="never") # noqa E1>
- ocean = cartopy.feature.OCEAN
-
- ax.add_feature(borders)
- ax.add_feature(coastlines)
- ax.add_feature(ocean, color="lightgrey")
+ ax.add_feature(cartopy.feature.OCEAN)
+ ax.add_feature(cartopy.feature.LAND)
+ ax.add_feature(cartopy.feature.COASTLINE)
+ ax.add_feature(cartopy.feature.BORDERS)
plt.tight_layout(pad=0)
@@ -208,7 +194,8 @@ def proj_area_attrs_section(area: 'geom.AreaDefinition') -> str: # noqa F821
"""
resolution_str = "/".join([str(round(x, 1)) for x in area.resolution])
- proj_dict = area.proj_dict
+ with ignore_pyproj_proj_warnings():
+ proj_dict = area.proj_dict
proj_str = "{{{}}}".format(", ".join(["'%s': '%s'" % (str(k), str(proj_dict[k])) for k in
sorted(proj_dict.keys())]))
area_units = proj_dict.get("units", "")
@@ -222,7 +209,7 @@ def proj_area_attrs_section(area: 'geom.AreaDefinition') -> str: # noqa F821
f"<dt>Width/Height</dt><dd>{area.width}/{area.height} Pixel</dd>"
f"<dt>Resolution x/y (SSP)</dt><dd>{resolution_str} {area_units}</dd>"
f"<dt>Extent (ll_x, ll_y, ur_x, ur_y)</dt>"
- f"<dd>{tuple(round(x, 4) for x in area.area_extent)}</dd>"
+ f"<dd>{tuple(round(float(x), 4) for x in area.area_extent)}</dd>"
"</dl>"
)
@@ -244,25 +231,30 @@ def swath_area_attrs_section(area: 'geom.SwathDefinition') -> str: # noqa F821
- Improve resolution estimation from lat/lon arrays. Maybe use CoordinateDefinition.geocentric_resolution?
"""
- if isinstance(area.lons, np.ndarray) & isinstance(area.lats, np.ndarray):
- # Calculate and estimated resolution from lats/lons in meter
- area_name = "Arbitrary Swath"
- resolution_y = np.mean(area.lats[0:-1, :] - area.lats[1::, :])
- resolution_x = np.mean(area.lons[:, 1::] - area.lons[:, 0:-1])
- resolution = np.mean(np.array([resolution_x, resolution_y]))
- resolution = np.round(40075000 * resolution / 360, 1)
- resolution_str = f"{resolution}/{resolution}"
+ if np.ndim(area.lons) == 1:
+ area_name = "1D Swath"
+ resolution_str = "NAxNA"
+ height, width = "NA", "NA"
area_units = "m"
else:
- lon_attrs = area.lons.attrs
- lat_attrs = area.lats.attrs
+ if isinstance(area.lons, np.ndarray) and isinstance(area.lats, np.ndarray):
+ # Calculate and estimated resolution from lats/lons in meter
+ area_name = "Arbitrary Swath"
+ resolution_y = np.mean(area.lats[0:-1, :] - area.lats[1::, :])
+ resolution_x = np.mean(area.lons[:, 1::] - area.lons[:, 0:-1])
+ resolution = np.mean(np.array([resolution_x, resolution_y]))
+ resolution = np.round(40075000 * resolution / 360, 1)
+ resolution_str = f"{resolution}x{resolution}"
+ else:
+ lon_attrs = area.lons.attrs
+ lat_attrs = area.lats.attrs
+
+ # use resolution from lat/lons dataarray attributes -> are these always set? -> Maybe try/except?
+ area_name = f"{lon_attrs.get('sensor')} Swath"
+ resolution_str = "x".join([str(round(x.get("resolution"), 1)) for x in [lat_attrs, lon_attrs]])
- # use resolution from lat/lons dataarray attributes -> are these always set? -> Maybe try/except?
- area_name = f"{lon_attrs.get('sensor')} swath"
- resolution_str = "/".join([str(round(x.get("resolution"), 1)) for x in [lat_attrs, lon_attrs]])
area_units = "m"
-
- height, width = area.lons.shape
+ height, width = area.lons.shape
attrs_icon = _icon("icon-file-text2")
@@ -270,7 +262,7 @@ def swath_area_attrs_section(area: 'geom.SwathDefinition') -> str: # noqa F821
# f"<dt>Area name</dt><dd>{area_name}</dd>"
f"<dt>Description</dt><dd>{area_name}</dd>"
f"<dt>Width/Height</dt><dd>{width}/{height} Pixel</dd>"
- f"<dt>Resolution x/y (SSP)</dt><dd>{resolution_str} {area_units}</dd>"
+ f"<dt>Resolution XxY (SSP)</dt><dd>{resolution_str} {area_units}</dd>"
"</dl>"
)
@@ -314,14 +306,15 @@ def swath_area_attrs_section(area: 'geom.SwathDefinition') -> str: # noqa F821
return f"{coll}"
-def area_repr(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'], include_header: Optional[bool] = True, # noqa F821
- include_static_files: Optional[bool] = True):
+def area_repr(area: Union['geom.AreaDefinition', 'geom.SwathDefinition'],
+ include_header: bool = True,
+ include_static_files: bool = True):
"""Return html repr of an AreaDefinition.
Args:
area : Area definition.
include_header : If true a header with object type will be included in
- the html. This is mainly intented for display in Jupyter Notebooks. For the
+ the html. This is mainly intended for display in Jupyter Notebooks. For the
display in the overview of area definitions for the Satpy documentation this
should be set to false.
include_static_files : Load and include css and html needed for representation.
=====================================
pyresample/bucket/__init__.py
=====================================
@@ -68,7 +68,7 @@ def _expand_bin_statistics(bins, unique_bin, unique_idx, weights_sorted):
# assign the valid index to array
weight_idx[unique_bin[~unique_bin.mask].data] = unique_idx[~unique_bin.mask]
- return weights_sorted[weight_idx] # last value of weigths_sorted always nan
+ return weights_sorted[weight_idx] # last value of weights_sorted always nan
@dask.delayed(pure=True)
@@ -202,20 +202,29 @@ class BucketResampler(object):
target_shape = self.target_area.shape
self.idxs = self.y_idxs * target_shape[1] + self.x_idxs
- def get_sum(self, data, skipna=True):
+ def get_sum(self, data, fill_value=np.nan, skipna=True, empty_bucket_value=0):
"""Calculate sums for each bin with drop-in-a-bucket resampling.
Parameters
----------
data : Numpy or Dask array
Data to be binned and summed.
+ fill_value : float
+ Fill value of the input data marking missing/invalid values.
+ Default: np.nan
skipna : boolean (optional)
- If True, skips NaN values for the sum calculation
- (similarly to Numpy's `nansum`). Buckets containing only NaN are set to zero.
- If False, sets the bucket to NaN if one or more NaN values are present in the bucket
- (similarly to Numpy's `sum`).
- In both cases, empty buckets are set to 0.
- Default: True
+ If True, skips missing values (as marked by NaN or `fill_value`) for the sum calculation
+ (similarly to Numpy's `nansum`). Buckets containing only missing values are set to `empty_bucket_value`.
+ If False, sets the bucket to fill_value if one or more missing values are present in the bucket
+ (similarly to Numpy's `sum`).
+ In both cases, empty buckets are set to `empty_bucket_value`.
+ Default: True
+ empty_bucket_value : float
+ Set empty buckets to the given value. Empty buckets are considered as the buckets with value 0.
+ Note that a bucket could become 0 as the result of a sum
+ of positive and negative values. If the user needs to identify these zero-buckets reliably,
+ `get_count()` can be used for this purpose.
+ Default: 0
Returns
-------
@@ -228,8 +237,9 @@ class BucketResampler(object):
data = data.data
data = data.ravel()
- # Remove NaN values from the data when used as weights
- weights = da.where(np.isnan(data), 0, data)
+ # Remove fill_values values from the data when used as weights
+ invalid_mask = _get_invalid_mask(data, fill_value)
+ weights = da.where(invalid_mask, 0, data)
# Rechunk indices to match the data chunking
if weights.chunks != self.idxs.chunks:
@@ -241,16 +251,19 @@ class BucketResampler(object):
weights=weights, density=False)
# TODO remove following line in favour of weights = data when dask histogram bug (issue #6935) is fixed
- sums = self._mask_bins_with_nan_if_not_skipna(skipna, data, out_size, sums)
+ sums = self._mask_bins_with_nan_if_not_skipna(skipna, data, out_size, sums, fill_value)
+
+ if empty_bucket_value != 0:
+ sums = da.where(sums == 0, empty_bucket_value, sums)
return sums.reshape(self.target_area.shape)
- def _mask_bins_with_nan_if_not_skipna(self, skipna, data, out_size, statistic):
+ def _mask_bins_with_nan_if_not_skipna(self, skipna, data, out_size, statistic, fill_value):
if not skipna:
- nans = np.isnan(data)
- nan_bins, _ = da.histogram(self.idxs[nans], bins=out_size,
- range=(0, out_size))
- statistic = da.where(nan_bins > 0, np.nan, statistic)
+ missing_val = _get_invalid_mask(data, fill_value)
+ missing_val_bins, _ = da.histogram(self.idxs[missing_val], bins=out_size,
+ range=(0, out_size))
+ statistic = da.where(missing_val_bins > 0, fill_value, statistic)
return statistic
def _call_bin_statistic(self, statistic_method, data, fill_value=None, skipna=None):
@@ -456,6 +469,14 @@ class BucketResampler(object):
return results
+def _get_invalid_mask(data, fill_value):
+ """Get a boolean array where values equal to fill_value in data are True."""
+ if np.isnan(fill_value):
+ return np.isnan(data)
+ else:
+ return data == fill_value
+
+
def round_to_resolution(arr, resolution):
"""Round the values in *arr* to closest resolution element.
=====================================
pyresample/geometry.py
=====================================
@@ -1451,11 +1451,14 @@ def masked_ints(func):
is_scalar = np.isscalar(xm) and np.isscalar(ym)
x__, y__ = func(self, xm, ym)
+ epsilon = 0.02 # arbitrary buffer for floating point precision
+ x_mask = ((x__ < -0.5 - epsilon) | (x__ > self.width - 0.5 + epsilon))
+ y_mask = ((y__ < -0.5 - epsilon) | (y__ > self.height - 0.5 + epsilon))
+ x__ = np.clip(x__, 0, self.width - 1)
+ y__ = np.clip(y__, 0, self.height - 1)
x__ = np.round(x__).astype(int)
y__ = np.round(y__).astype(int)
- x_mask = ((x__ < 0) | (x__ >= self.width))
- y_mask = ((y__ < 0) | (y__ >= self.height))
x_masked = np.ma.masked_array(x__, mask=x_mask, copy=False)
y_masked = np.ma.masked_array(y__, mask=y_mask, copy=False)
if is_scalar:
=====================================
pyresample/static/css/style.css
=====================================
@@ -29,7 +29,7 @@ body.vscode-dark {
.pyresample-wrap {
display: block !important;
min-width: 700px;
- max-width: 900px;
+ max-width: 1050px;
}
.pyresample-text-repr-fallback {
@@ -90,19 +90,6 @@ body.vscode-dark {
margin-right: 0.5em;
}
-/* und noch der Doppelpunkt */
-.pyresample-area-section-details dt:after {
- content: ": ";
-}
-
-/* ein Clearfix für Folge-dd-Elemente */
-.pyresample-area-section-details dd:after {
- clear: left;
- content: " ";
- display: block;
-}
-/* end attribute list formatting */
-
.pyresample-area-sections {
display: grid;
grid-template-columns: auto;
=====================================
pyresample/test/test_bucket.py
=====================================
@@ -17,328 +17,377 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Test the bucket resampler."""
-import unittest
from unittest.mock import MagicMock, patch
import dask
import dask.array as da
import numpy as np
+import pytest
import xarray as xr
from pyresample import bucket, create_area_def
+from pyresample.bucket import _get_invalid_mask
from pyresample.geometry import AreaDefinition
from pyresample.test.utils import CustomScheduler
+CHUNKS = 2
+WIDTH = 2560
+HEIGHT = 2048
-class Test(unittest.TestCase):
- """Test bucket resampler."""
- adef = AreaDefinition('eurol', 'description', '',
+ at pytest.fixture(scope="module")
+def adef():
+ """Get AreaDefinition for tests."""
+ return AreaDefinition('eurol',
+ 'description',
+ '',
{'ellps': 'WGS84',
'lat_0': '90.0',
'lat_ts': '60.0',
'lon_0': '0.0',
- 'proj': 'stere'}, 2560, 2048,
+ 'proj': 'stere'},
+ 2560,
+ 2048,
(-3780000.0, -7644000.0, 3900000.0, -1500000.0))
- chunks = 2
- lons = da.from_array(np.array([[25., 25.], [25., 25.]]),
- chunks=chunks)
- lats = da.from_array(np.array([[60., 60.00001], [60.2, 60.3]]),
- chunks=chunks)
-
- def setUp(self):
- self.resampler = bucket.BucketResampler(self.adef, self.lons, self.lats)
-
- @patch('pyresample.bucket.Proj')
- @patch('pyresample.bucket.BucketResampler._get_indices')
- def test_init(self, get_indices, prj):
- resampler = bucket.BucketResampler(self.adef, self.lons, self.lats)
- get_indices.assert_called_once()
- prj.assert_called_once_with(self.adef.proj_dict)
- self.assertTrue(hasattr(resampler, 'target_area'))
- self.assertTrue(hasattr(resampler, 'source_lons'))
- self.assertTrue(hasattr(resampler, 'source_lats'))
- self.assertTrue(hasattr(resampler, 'x_idxs'))
- self.assertTrue(hasattr(resampler, 'y_idxs'))
- self.assertTrue(hasattr(resampler, 'idxs'))
- self.assertTrue(hasattr(resampler, 'get_sum'))
- self.assertTrue(hasattr(resampler, 'get_count'))
- self.assertTrue(hasattr(resampler, 'get_min'))
- self.assertTrue(hasattr(resampler, 'get_max'))
- self.assertTrue(hasattr(resampler, 'get_abs_max'))
- self.assertTrue(hasattr(resampler, 'get_average'))
- self.assertTrue(hasattr(resampler, 'get_fractions'))
- self.assertIsNone(resampler.counts)
-
- def test_round_to_resolution(self):
- """Test rounding to given resolution."""
- # Scalar, integer resolution
- self.assertEqual(bucket.round_to_resolution(5.5, 2.), 6)
- # Scalar, non-integer resolution
- self.assertEqual(bucket.round_to_resolution(5.5, 1.7), 5.1)
- # List
- self.assertTrue(np.all(bucket.round_to_resolution([4.2, 5.6], 2) ==
- np.array([4., 6.])))
- # Numpy array
- self.assertTrue(np.all(bucket.round_to_resolution(np.array([4.2, 5.6]), 2) ==
- np.array([4., 6.])))
- # Dask array
- self.assertTrue(
- np.all(bucket.round_to_resolution(da.array([4.2, 5.6]), 2) ==
- np.array([4., 6.])))
-
- def test_get_proj_coordinates(self):
- """Test calculation of projection coordinates."""
- prj = MagicMock()
- prj.return_value = ([3.1, 3.1, 3.1], [4.8, 4.8, 4.8])
- lons = [1., 1., 1.]
- lats = [2., 2., 2.]
- self.resampler.prj = prj
- result = self.resampler._get_proj_coordinates(lons, lats)
- prj.assert_called_once_with(lons, lats)
- self.assertTrue(isinstance(result, np.ndarray))
- np.testing.assert_equal(result, np.array([[3.1, 3.1, 3.1],
- [4.8, 4.8, 4.8]]))
-
- def test_get_bucket_indices(self):
- """Test calculation of array indices."""
- # Ensure nothing is calculated
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- self.resampler._get_indices()
- x_idxs, y_idxs = da.compute(self.resampler.x_idxs,
- self.resampler.y_idxs)
- np.testing.assert_equal(x_idxs, np.array([1710, 1710, 1707, 1705]))
- np.testing.assert_equal(y_idxs, np.array([465, 465, 459, 455]))
-
- # Additional small test case
- adef = create_area_def(
- area_id='test',
- projection={'proj': 'latlong'},
- width=2, height=2,
- center=(0, 0),
- resolution=10)
- lons = da.from_array(
- np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, -10.1, 0]),
- chunks=2)
- lats = da.from_array(
- np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, 0, 10.1]),
- chunks=2)
- resampler = bucket.BucketResampler(source_lats=lats,
- source_lons=lons,
- target_area=adef)
+
+
+ at pytest.fixture(scope="module")
+def lons():
+ """Get longitudes for tests."""
+ return da.from_array(np.array([[25., 25.], [25., 25.]]), chunks=CHUNKS)
+
+
+ at pytest.fixture(scope="module")
+def lats():
+ """Get latitudes for tests."""
+ return da.from_array(np.array([[60., 60.00001], [60.2, 60.3]]), chunks=CHUNKS)
+
+
+ at pytest.fixture(scope="module")
+def resampler(adef, lons, lats):
+ """Get initialised resampler for tests."""
+ return bucket.BucketResampler(adef, lons, lats)
+
+
+ at patch('pyresample.bucket.Proj')
+ at patch('pyresample.bucket.BucketResampler._get_indices')
+def test_init(get_indices, prj, adef, lons, lats):
+ """Test the init method of the BucketResampler."""
+ resampler = bucket.BucketResampler(adef, lons, lats)
+
+ get_indices.assert_called_once()
+ prj.assert_called_once_with(adef.proj_dict)
+
+ assert hasattr(resampler, 'target_area')
+ assert hasattr(resampler, 'source_lons')
+ assert hasattr(resampler, 'source_lats')
+ assert hasattr(resampler, 'x_idxs')
+ assert hasattr(resampler, 'y_idxs')
+ assert hasattr(resampler, 'idxs')
+ assert hasattr(resampler, 'get_sum')
+ assert hasattr(resampler, 'get_count')
+ assert hasattr(resampler, 'get_min')
+ assert hasattr(resampler, 'get_max')
+ assert hasattr(resampler, 'get_abs_max')
+ assert hasattr(resampler, 'get_average')
+ assert hasattr(resampler, 'get_fractions')
+ assert resampler.counts is None
+
+
+def test_round_to_resolution():
+ """Test rounding to given resolution."""
+ # Scalar, integer resolution
+ assert bucket.round_to_resolution(5.5, 2.) == 6
+ # Scalar, non-integer resolution
+ assert bucket.round_to_resolution(5.5, 1.7) == 5.1
+ # List
+ assert np.all(bucket.round_to_resolution([4.2, 5.6], 2) == np.array([4., 6.]))
+ # Numpy array
+ assert np.all(bucket.round_to_resolution(np.array([4.2, 5.6]), 2) == np.array([4., 6.]))
+ # Dask array
+ assert np.all(bucket.round_to_resolution(da.array([4.2, 5.6]), 2) == np.array([4., 6.]))
+
+
+def test_get_proj_coordinates(adef, lons, lats):
+ """Test calculation of projection coordinates."""
+ resampler = bucket.BucketResampler(source_lats=lats, source_lons=lons, target_area=adef)
+ prj = MagicMock()
+ prj.return_value = ([3.1, 3.1, 3.1], [4.8, 4.8, 4.8])
+ lons = [1., 1., 1.]
+ lats = [2., 2., 2.]
+ resampler.prj = prj
+
+ result = resampler._get_proj_coordinates(lons, lats)
+
+ prj.assert_called_once_with(lons, lats)
+ assert isinstance(result, np.ndarray)
+ np.testing.assert_equal(result, np.array([[3.1, 3.1, 3.1],
+ [4.8, 4.8, 4.8]]))
+
+
+def test_get_bucket_indices(resampler):
+ """Test calculation of array indices."""
+ # Ensure nothing is calculated
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
resampler._get_indices()
- np.testing.assert_equal(resampler.x_idxs, np.array([-1, 0, 0, 1, 1, 1, -1, -1, -1]))
- np.testing.assert_equal(resampler.y_idxs, np.array([-1, 1, 1, 1, 0, 0, -1, -1, -1]))
-
- def _get_sum_result(self, data, **kwargs):
- """Compute the bucket average with kwargs and check that no dask computation is performed."""
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_sum(data, **kwargs)
- return result.compute()
-
- def test_get_sum_valid_data(self):
- """Test drop-in-a-bucket sum for valid data input."""
- data = da.from_array(np.array([[2., 3.], [7., 16.]]),
- chunks=self.chunks)
-
- result = self._get_sum_result(data)
-
- # first two values are in same bin
- self.assertEqual(np.count_nonzero(result == 5), 1)
- # others are in separate bins
- self.assertEqual(np.count_nonzero(result == 7), 1)
- self.assertEqual(np.count_nonzero(result == 16), 1)
-
- self.assertEqual(result.shape, self.adef.shape)
-
- # Test that also xarray.DataArrays work (same output)
- data = xr.DataArray(data)
- np.testing.assert_array_equal(result, self._get_sum_result(data))
-
- def test_get_sum_nan_data_skipna_false(self):
- """Test drop-in-a-bucket sum for data input with nan and skipna False."""
- data = da.from_array(np.array([[2., np.nan], [5., np.nan]]),
- chunks=self.chunks)
-
- result = self._get_sum_result(data, skipna=False)
- # 2 + nan is nan, all-nan bin is nan
- self.assertEqual(np.count_nonzero(np.isnan(result)), 2)
- # rest is 0
- self.assertEqual(np.nanmin(result), 0)
-
- def test_get_sum_nan_data_skipna_true(self):
- """Test drop-in-a-bucket sum for data input with nan and skipna True."""
- data = da.from_array(np.array([[2., np.nan], [5., np.nan]]),
- chunks=self.chunks)
-
- result = self._get_sum_result(data, skipna=True)
- # 2 + nan is 2
- self.assertEqual(np.count_nonzero(result == 2.), 1)
- # all-nan and rest is 0
- self.assertEqual(np.count_nonzero(np.isnan(result)), 0)
- self.assertEqual(np.nanmin(result), 0)
-
- def test_get_count(self):
- """Test drop-in-a-bucket sum."""
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_count()
- result = result.compute()
- self.assertTrue(np.max(result) == 2)
- self.assertEqual(np.sum(result == 1), 2)
- self.assertEqual(np.sum(result == 2), 1)
- self.assertTrue(self.resampler.counts is not None)
-
- def _get_min_result(self, data, **kwargs):
- """Compute the bucket average with kwargs and check that no dask computation is performed."""
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_min(data, **kwargs)
- return result.compute()
-
- def test_get_min(self):
- """Test min bucket resampling."""
- data = da.from_array(np.array([[2, 11], [5, np.nan]]),
- chunks=self.chunks)
- result = self._get_min_result(data)
- # test multiple entries minimum
- self.assertEqual(np.count_nonzero(result == 2), 1)
- # test single entry minimum
- self.assertEqual(np.count_nonzero(result == 5), 1)
- # test that minimum of bucket with only nan is nan, and empty buckets are nan
- self.assertEqual(np.count_nonzero(~np.isnan(result)), 2)
-
- def _get_max_result(self, data, **kwargs):
- """Compute the bucket max with kwargs and check that no dask computation is performed."""
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_max(data, **kwargs)
- return result.compute()
-
- def test_get_max(self):
- """Test max bucket resampling."""
- data = da.from_array(np.array([[2, 11], [5, np.nan]]),
- chunks=self.chunks)
- result = self._get_max_result(data)
- # test multiple entries maximum
- self.assertEqual(np.count_nonzero(result == 11), 1)
- # test single entry maximum
- self.assertEqual(np.count_nonzero(result == 5), 1)
- # test that minimum of bucket with only nan is nan, and empty buckets are nan
- self.assertEqual(np.count_nonzero(~np.isnan(result)), 2)
-
- def test_get_abs_max(self):
- """Test abs max bucket resampling."""
- data = da.from_array(np.array([[2, -11], [5, np.nan]]),
- chunks=self.chunks)
- result = self._get_abs_max_result(data)
- # test multiple entries absolute maximum
- self.assertEqual(np.count_nonzero(result == -11), 1)
- # test single entry maximum
- self.assertEqual(np.count_nonzero(result == 5), 1)
- # test that minimum of bucket with only nan is nan, and empty buckets are nan
- self.assertEqual(np.count_nonzero(~np.isnan(result)), 2)
-
- def _get_abs_max_result(self, data, **kwargs):
- """Compute the bucket abs max with kwargs and check that no dask computation is performed."""
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_abs_max(data, **kwargs)
- return result.compute()
-
- def _get_average_result(self, data, **kwargs):
- """Compute the bucket average with kwargs and check that no dask computation is performed."""
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_average(data, **kwargs)
- return result.compute()
-
- def test_get_average_basic(self):
- """Test averaging bucket resampling."""
- data = da.from_array(np.array([[2, 11], [5, np.nan]]),
- chunks=self.chunks)
- result = self._get_average_result(data)
- # test multiple entries average
- self.assertEqual(np.count_nonzero(result == 6.5), 1)
- # test single entry average
- self.assertEqual(np.count_nonzero(result == 5), 1)
- # test that average of bucket with only nan is nan, and empty buckets are nan
- self.assertEqual(np.count_nonzero(~np.isnan(result)), 2)
-
- def test_get_average_with_fillvalue_for_output(self):
- """Test averaging bucket resampling with defined fill_value for output."""
- data = da.from_array(np.array([[2, 11], [5, np.nan]]),
- chunks=self.chunks)
- # test fill_value other than np.nan
- result = self._get_average_result(data, fill_value=-1)
- # check that all empty buckets are fill_value
- self.assertEqual(np.count_nonzero(result != -1), 2)
-
- def test_get_average_skipna_true(self):
- """Test averaging bucket resampling with skipna True."""
- # test skipna
- data = da.from_array(np.array([[2, np.nan], [np.nan, np.nan]]),
- chunks=self.chunks)
- result = self._get_average_result(data, skipna=True)
- # test that average of 2 and np.nan is 2 for skipna=True
- self.assertEqual(np.count_nonzero(result == 2), 1)
-
- def test_get_average_skipna_false(self):
- """Test averaging bucket resampling with skipna False."""
- data = da.from_array(np.array([[2, np.nan], [np.nan, np.nan]]),
- chunks=self.chunks)
- result = self._get_average_result(data, skipna=False)
- # test that average of 2 and np.nan is nan for skipna=False
- self.assertTrue(np.all(np.isnan(result)))
-
- def test_get_average_only_nan_input(self):
- """Test averaging bucket resampling with only NaN as input."""
- data = da.from_array(np.array([[np.nan, np.nan], [np.nan, np.nan]]),
- chunks=self.chunks)
- result = self._get_average_result(data, skipna=True)
- # test that average of np.nan and np.nan is np.nan for both skipna
- self.assertTrue(np.all(np.isnan(result)))
- np.testing.assert_array_equal(result, self._get_average_result(data, skipna=False))
-
- def test_get_average_with_fill_value_in_input(self):
- """Test averaging bucket resampling with fill_value in input and skipna True."""
- # test that fill_value in input is recognised as missing value
- data = da.from_array(np.array([[2, -1], [-1, np.nan]]),
- chunks=self.chunks)
- result = self._get_average_result(data, fill_value=-1, skipna=True)
- # test that average of 2 and -1 (missing value) is 2
- self.assertEqual(np.count_nonzero(result == 2), 1)
- # test than all other buckets are -1
- self.assertEqual(np.count_nonzero(result != -1), 1)
-
- def test_resample_bucket_fractions(self):
- """Test fraction calculations for categorical data."""
- data = da.from_array(np.array([[2, 4], [2, 2]]),
- chunks=self.chunks)
- categories = [1, 2, 3, 4]
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_fractions(data, categories=categories)
- self.assertEqual(set(categories), set(result.keys()))
- res = result[1].compute()
- self.assertTrue(np.nanmax(res) == 0.)
- res = result[2].compute()
- self.assertTrue(np.nanmax(res) == 1.)
- self.assertTrue(np.nanmin(res) == 0.5)
- res = result[3].compute()
- self.assertTrue(np.nanmax(res) == 0.)
- res = result[4].compute()
- self.assertTrue(np.nanmax(res) == 0.5)
- self.assertTrue(np.nanmin(res) == 0.)
- # There should be NaN values
- self.assertTrue(np.any(np.isnan(res)))
-
- # Use a fill value
- with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
- result = self.resampler.get_fractions(data, categories=categories,
- fill_value=-1)
-
- # There should not be any NaN values
- for i in categories:
- res = result[i].compute()
- self.assertFalse(np.any(np.isnan(res)))
- self.assertTrue(np.min(res) == -1)
-
- # No categories given, need to compute the data once to get
- # the categories
- with dask.config.set(scheduler=CustomScheduler(max_computes=1)):
- _ = self.resampler.get_fractions(data, categories=None)
+ x_idxs, y_idxs = da.compute(resampler.x_idxs, resampler.y_idxs)
+ np.testing.assert_equal(x_idxs, np.array([1710, 1710, 1707, 1705]))
+ np.testing.assert_equal(y_idxs, np.array([465, 465, 459, 455]))
+
+
+def test_get_bucket_indices_on_latlong():
+ """Test calculation of array indices on latlong grid."""
+ adef = create_area_def(
+ area_id='test',
+ projection={'proj': 'latlong'},
+ width=2, height=2,
+ center=(0, 0),
+ resolution=10)
+ lons = da.from_array(np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, -10.1, 0]), chunks=CHUNKS)
+ lats = da.from_array(np.array([-10.0, -9.9, -0.1, 0, 0.1, 9.9, 10.0, 0, 10.1]), chunks=CHUNKS)
+ resampler = bucket.BucketResampler(source_lats=lats, source_lons=lons, target_area=adef)
+ resampler._get_indices()
+
+ np.testing.assert_equal(resampler.x_idxs, np.array([-1, 0, 0, 1, 1, 1, -1, -1, -1]))
+ np.testing.assert_equal(resampler.y_idxs, np.array([-1, 1, 1, 1, 0, 0, -1, -1, -1]))
+
+
+def _get_sum_result(resampler, data, **kwargs):
+ """Compute the bucket average with kwargs and check that no dask computation is performed."""
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_sum(data, **kwargs)
+ return result.compute()
+
+
+def test_get_sum_valid_data(resampler, adef):
+ """Test drop-in-a-bucket sum for valid data input."""
+ data = da.from_array(np.array([[2., 3.], [7., 16.]]), chunks=CHUNKS)
+ result = _get_sum_result(resampler, data)
+
+ # first two values are in same bin
+ assert np.count_nonzero(result == 5) == 1
+ # others are in separate bins
+ assert np.count_nonzero(result == 7) == 1
+ assert np.count_nonzero(result == 16) == 1
+ assert result.shape == adef.shape
+
+ # Test that also xarray.DataArrays work (same output)
+ data = xr.DataArray(data)
+ np.testing.assert_array_equal(result, _get_sum_result(resampler, data))
+
+
+def _equal_or_both_nan(val1, val2):
+ return val1 == val2 or (np.isnan(val1) and np.isnan(val2))
+
+
+ at pytest.mark.parametrize("skipna", [True, False])
+ at pytest.mark.parametrize("fill_value", [np.nan, 255, -1])
+ at pytest.mark.parametrize("empty_bucket_value", [0, 4095, np.nan, -1])
+def test_get_sum_skipna_fillvalue_empty_bucket_value(resampler, skipna, fill_value, empty_bucket_value):
+ """Test drop-in-a-bucket sum for invalid data input and according arguments."""
+ data = da.from_array(np.array([[2., fill_value], [5., fill_value]]), chunks=CHUNKS)
+ result = _get_sum_result(resampler, data,
+ skipna=skipna,
+ fill_value=fill_value,
+ empty_bucket_value=empty_bucket_value)
+ n_target_bkt = WIDTH * HEIGHT
+
+ # 5 is untouched in a single bin, in any case
+ n_bkt_with_val_5 = 1
+
+ if skipna:
+ # 2 + fill_value is 2 (nansum)
+ n_bkt_with_val_2 = 1
+ # and fill_value+fill_value is empty_bucket_value,
+ # hence no fill_value bkt are left
+ n_bkt_with_val_fill_value = 0
+ else:
+ # 2 + fill_value is fill_value (sum)
+ n_bkt_with_val_2 = 0
+ # and fill_value + fill_value is fill_value, so
+ n_bkt_with_val_fill_value = 2
+
+ n_bkt_with_empty_value = n_target_bkt - n_bkt_with_val_fill_value - n_bkt_with_val_5 - n_bkt_with_val_2
+
+ # special case
+ if _equal_or_both_nan(fill_value, empty_bucket_value):
+ # the fill and empty values are equal, so they should be added up
+ n_bkt_with_empty_value = n_bkt_with_val_fill_value = n_bkt_with_empty_value + n_bkt_with_val_fill_value
+
+ assert np.count_nonzero(result == 5.) == n_bkt_with_val_5
+ assert np.count_nonzero(result == 2.) == n_bkt_with_val_2
+ assert np.count_nonzero(_get_invalid_mask(result, fill_value)) == n_bkt_with_val_fill_value
+ assert np.count_nonzero(_get_invalid_mask(result, empty_bucket_value)) == n_bkt_with_empty_value
+
+
+def test_get_count(resampler):
+ """Test drop-in-a-bucket sum."""
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_count()
+ result = result.compute()
+ assert np.max(result) == 2
+ assert np.sum(result == 1) == 2
+ assert np.sum(result == 2) == 1
+ assert resampler.counts is not None
+
+
+def _get_min_result(resampler, data, **kwargs):
+ """Compute the bucket average with kwargs and check that no dask computation is performed."""
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_min(data, **kwargs)
+ return result.compute()
+
+
+def test_get_min(resampler):
+ """Test min bucket resampling."""
+ data = da.from_array(np.array([[2, 11], [5, np.nan]]), chunks=CHUNKS)
+ result = _get_min_result(resampler, data)
+ # test multiple entries minimum
+ assert np.count_nonzero(result == 2) == 1
+ # test single entry minimum
+ assert np.count_nonzero(result == 5) == 1
+ # test that minimum of bucket with only nan is nan, and empty buckets are nan
+ assert np.count_nonzero(~np.isnan(result)) == 2
+
+
+def _get_max_result(resampler, data, **kwargs):
+ """Compute the bucket max with kwargs and check that no dask computation is performed."""
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_max(data, **kwargs)
+ return result.compute()
+
+
+def test_get_max(resampler):
+ """Test max bucket resampling."""
+ data = da.from_array(np.array([[2, 11], [5, np.nan]]), chunks=CHUNKS)
+ result = _get_max_result(resampler, data)
+ # test multiple entries maximum
+ assert np.count_nonzero(result == 11) == 1
+ # test single entry maximum
+ assert np.count_nonzero(result == 5) == 1
+ # test that minimum of bucket with only nan is nan, and empty buckets are nan
+ assert np.count_nonzero(~np.isnan(result)) == 2
+
+
+def _get_abs_max_result(resampler, data, **kwargs):
+ """Compute the bucket abs max with kwargs and check that no dask computation is performed."""
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_abs_max(data, **kwargs)
+ return result.compute()
+
+
+def test_get_abs_max(resampler):
+ """Test abs max bucket resampling."""
+ data = da.from_array(np.array([[2, -11], [5, np.nan]]), chunks=CHUNKS)
+ result = _get_abs_max_result(resampler, data)
+ # test multiple entries absolute maximum
+ assert np.count_nonzero(result == -11) == 1
+ # test single entry maximum
+ assert np.count_nonzero(result == 5) == 1
+ # test that minimum of bucket with only nan is nan, and empty buckets are nan
+ assert np.count_nonzero(~np.isnan(result)) == 2
+
+
+def _get_average_result(resampler, data, **kwargs):
+ """Compute the bucket average with kwargs and check that no dask computation is performed."""
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_average(data, **kwargs)
+ return result.compute()
+
+
+def test_get_average_basic(resampler):
+ """Test averaging bucket resampling."""
+ data = da.from_array(np.array([[2, 11], [5, np.nan]]), chunks=CHUNKS)
+ result = _get_average_result(resampler, data)
+ # test multiple entries average
+ assert np.count_nonzero(result == 6.5) == 1
+ # test single entry average
+ assert np.count_nonzero(result == 5) == 1
+ # test that average of bucket with only nan is nan, and empty buckets are nan
+ assert np.count_nonzero(~np.isnan(result)) == 2
+
+
+def test_get_average_with_fillvalue_for_output(resampler):
+ """Test averaging bucket resampling with defined fill_value for output."""
+ data = da.from_array(np.array([[2, 11], [5, np.nan]]), chunks=CHUNKS)
+ # test fill_value other than np.nan
+ result = _get_average_result(resampler, data, fill_value=-1)
+ # check that all empty buckets are fill_value
+ assert np.count_nonzero(result != -1) == 2
+
+
+def test_get_average_skipna_true(resampler):
+ """Test averaging bucket resampling with skipna True."""
+ # test skipna
+ data = da.from_array(np.array([[2, np.nan], [np.nan, np.nan]]), chunks=CHUNKS)
+ result = _get_average_result(resampler, data, skipna=True)
+ # test that average of 2 and np.nan is 2 for skipna=True
+ assert np.count_nonzero(result == 2) == 1
+
+
+def test_get_average_skipna_false(resampler):
+ """Test averaging bucket resampling with skipna False."""
+ data = da.from_array(np.array([[2, np.nan], [np.nan, np.nan]]), chunks=CHUNKS)
+ result = _get_average_result(resampler, data, skipna=False)
+ # test that average of 2 and np.nan is nan for skipna=False
+ assert np.all(np.isnan(result))
+
+
+def test_get_average_only_nan_input(resampler):
+ """Test averaging bucket resampling with only NaN as input."""
+ data = da.from_array(np.array([[np.nan, np.nan], [np.nan, np.nan]]), chunks=CHUNKS)
+ result = _get_average_result(resampler, data, skipna=True)
+ # test that average of np.nan and np.nan is np.nan for both skipna
+ assert np.all(np.isnan(result))
+ np.testing.assert_array_equal(result, _get_average_result(resampler, data, skipna=False))
+
+
+def test_get_average_with_fill_value_in_input(resampler):
+ """Test averaging bucket resampling with fill_value in input and skipna True."""
+ # test that fill_value in input is recognised as missing value
+ data = da.from_array(np.array([[2, -1], [-1, np.nan]]), chunks=CHUNKS)
+ result = _get_average_result(resampler, data, fill_value=-1, skipna=True)
+ # test that average of 2 and -1 (missing value) is 2
+ assert np.count_nonzero(result == 2) == 1
+ # test that all other buckets are -1
+ assert np.count_nonzero(result != -1) == 1
+
+
+def test_resample_bucket_fractions(resampler):
+ """Test fraction calculations for categorical data."""
+ data = da.from_array(np.array([[2, 4], [2, 2]]), chunks=CHUNKS)
+ categories = [1, 2, 3, 4]
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_fractions(data, categories=categories)
+ assert set(categories) == set(result.keys())
+
+ res = result[1].compute()
+ assert np.nanmax(res) == 0.
+
+ res = result[2].compute()
+ assert np.nanmax(res) == 1.
+ assert np.nanmin(res) == 0.5
+
+ res = result[3].compute()
+ assert np.nanmax(res) == 0.
+
+ res = result[4].compute()
+ assert np.nanmax(res) == 0.5
+ assert np.nanmin(res) == 0.
+ # There should be NaN values
+ assert np.any(np.isnan(res))
+
+ # Use a fill value
+ with dask.config.set(scheduler=CustomScheduler(max_computes=0)):
+ result = resampler.get_fractions(data, categories=categories, fill_value=-1)
+
+ # There should not be any NaN values
+ for i in categories:
+ res = result[i].compute()
+ assert not np.any(np.isnan(res))
+ assert np.min(res) == -1
+
+ # No categories given, need to compute the data once to get
+ # the categories
+ with dask.config.set(scheduler=CustomScheduler(max_computes=1)):
+ _ = resampler.get_fractions(data, categories=None)
=====================================
pyresample/test/test_geometry/test_area.py
=====================================
@@ -598,47 +598,6 @@ class TestAreaDefinition:
assert lons[0, 0] == -179.5
assert lats[0, 0] == 89.5
- def test_get_array_indices_from_lonlat_mask_actual_values(self, create_test_area):
- """Test that the masked values of get_array_indices_from_lonlat can be valid."""
- x_size = 3712
- y_size = 3712
- area_extent = (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)
- proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'lat_1': 25.,
- 'lat_2': 25., 'lon_0': 0.0, 'proj': 'lcc', 'units': 'm'}
- area_def = create_test_area(proj_dict, x_size, y_size, area_extent)
-
- # Choose a point just outside the area
- x, y = area_def.get_array_indices_from_lonlat([33.5], [-40.5])
- assert x.item() == 3723
- assert y.item() == 3746
-
- def test_lonlat2colrow(self, create_test_area):
- """Test lonlat2colrow."""
- x_size = 3712
- y_size = 3712
- area_extent = [-5570248.477339261, -5567248.074173444, 5567248.074173444, 5570248.477339261]
- proj_dict = {'a': '6378169.00',
- 'b': '6356583.80',
- 'h': '35785831.0',
- 'lon_0': '0.0',
- 'proj': 'geos'}
- area = create_test_area(proj_dict, x_size, y_size, area_extent)
-
- # Imatra, Wiesbaden
- longitudes = np.array([28.75242, 8.24932])
- latitudes = np.array([61.17185, 50.08258])
- cols__, rows__ = area.get_array_indices_from_lonlat(longitudes, latitudes)
-
- # test arrays
- cols_expects = np.array([2304, 2040])
- rows_expects = np.array([186, 341])
- np.testing.assert_array_equal(cols__, cols_expects)
- np.testing.assert_array_equal(rows__, rows_expects)
-
- # test scalars
- lon, lat = (-8.125547604568746, -14.345524111874646)
- assert area.get_array_indices_from_lonlat(lon, lat) == (1567, 2375)
-
def test_colrow2lonlat(self, create_test_area):
"""Test colrow2lonlat."""
# test square, symmetric areadef
@@ -800,8 +759,69 @@ class TestAreaDefinition:
actual = [(np.rad2deg(coord.lon), np.rad2deg(coord.lat)) for coord in area_def.corners]
np.testing.assert_allclose(actual, expected)
+ @pytest.mark.parametrize(
+ ("lon", "lat", "exp_mask", "exp_col", "exp_row"),
+ [
+ # Choose a point outside the area
+ (33.5, -40.5, True, 0.0, 0.0),
+ # A point just barely outside the left extent (near floating point precision)
+ (-63.62135, 37.253807, False, 0, 5),
+ # A point just barely outside the right extent (near floating point precision)
+ (63.59189, 37.26574, False, 3711, 5),
+ ]
+ )
+ def test_get_array_indices_from_lonlat_mask_actual_values(
+ self, create_test_area, lon, lat, exp_mask, exp_col, exp_row):
+ """Test masking behavior of get_array_indices_from_lonlat edge cases."""
+ x_size = 3712
+ y_size = 3712
+ area_extent = (-5570248.477339745, -5561247.267842293, 5567248.074173927, 5570248.477339745)
+ proj_dict = {'a': 6378169.0, 'b': 6356583.8, 'lat_1': 25.,
+ 'lat_2': 25., 'lon_0': 0.0, 'proj': 'lcc', 'units': 'm'}
+ area_def = create_test_area(proj_dict, x_size, y_size, area_extent)
+
+ x, y = area_def.get_array_indices_from_lonlat([lon], [lat])
+ if exp_mask:
+ assert x.mask.all()
+ assert y.mask.all()
+ else:
+ assert not x.mask.any()
+ assert not y.mask.any()
+ assert x.item() == exp_col
+ assert y.item() == exp_row
+
+ @pytest.mark.parametrize(
+ ("lons", "lats", "exp_cols", "exp_rows"),
+ [
+ # Imatra, Wiesbaden
+ (np.array([28.75242, 8.24932]), np.array([61.17185, 50.08258]),
+ np.array([2304, 2040]), np.array([186, 341])),
+ (-8.125547604568746, -14.345524111874646,
+ 1567, 2375),
+ ]
+ )
+ def test_get_array_indices_from_lonlat_geos(self, create_test_area, lons, lats, exp_cols, exp_rows):
+ """Test get_array_indices_from_lonlat."""
+ x_size = 3712
+ y_size = 3712
+ area_extent = [-5570248.477339261, -5567248.074173444, 5567248.074173444, 5570248.477339261]
+ proj_dict = {'a': '6378169.00',
+ 'b': '6356583.80',
+ 'h': '35785831.0',
+ 'lon_0': '0.0',
+ 'proj': 'geos'}
+ area = create_test_area(proj_dict, x_size, y_size, area_extent)
+
+ cols__, rows__ = area.get_array_indices_from_lonlat(lons, lats)
+ if hasattr(exp_cols, "shape"):
+ np.testing.assert_array_equal(cols__, exp_cols)
+ np.testing.assert_array_equal(rows__, exp_rows)
+ else:
+ assert cols__ == exp_cols
+ assert rows__ == exp_rows
+
def test_get_array_indices_from_lonlat(self, create_test_area):
- """Test the function get_xy_from_lonlat."""
+ """Test the function get_array_indices_from_lonlat."""
x_size = 2
y_size = 2
area_extent = [1000000, 0, 1050000, 50000]
=====================================
pyresample/utils/cartopy.py
=====================================
@@ -31,6 +31,7 @@ except ImportError:
DataArray = np.ndarray
import cartopy.crs as ccrs
+import shapely.geometry as sgeom
class Projection(ccrs.Projection):
@@ -69,3 +70,77 @@ class Projection(ccrs.Projection):
if self.bounds is not None:
x0, x1, y0, y1 = self.bounds
self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.
+
+ # For geostationary full disc the boundary needs to be the actuall boundary of the earth disc otherwise
+ # when ocean or land features are added to the cartopy plot these spill over.
+ if "Geostationary Satellite" not in crs.to_wkt():
+ self._boundary = super().boundary
+ else:
+ self._boundary = self._geos_boundary()
+
+ @staticmethod
+ def _geos_boundary():
+ """Calculate full disk boundary.
+
+ This code is copied over from the 'Geostationary' class in 'cartopy/lib/cartopy/crs.py'.
+ """
+ satellite_height = 35785831
+ false_easting = 0
+ false_northing = 0
+ a = float(ccrs.WGS84_SEMIMAJOR_AXIS)
+ b = float(ccrs.WGS84_SEMIMINOR_AXIS)
+ h = float(satellite_height)
+
+ # To find the bound we trace around where the line from the satellite
+ # is tangent to the surface. This involves trigonometry on a sphere
+ # centered at the satellite. The two scanning angles form two legs of
+ # triangle on this sphere--the hypotenuse "c" (angle arc) is controlled
+ # by distance from center to the edge of the ellipse being seen.
+
+ # This is one of the angles in the spherical triangle and used to
+ # rotate around and "scan" the boundary
+ angleA = np.linspace(0, -2 * np.pi, 91) # Clockwise boundary.
+
+ # Convert the angle around center to the proper value to use in the
+ # parametric form of an ellipse
+ th = np.arctan(a / b * np.tan(angleA))
+
+ # Given the position on the ellipse, what is the distance from center
+ # to the ellipse--and thus the tangent point
+ r = np.hypot(a * np.cos(th), b * np.sin(th))
+ sat_dist = a + h
+
+ # Using this distance, solve for sin and tan of c in the triangle that
+ # includes the satellite, Earth center, and tangent point--we need to
+ # figure out the location of this tangent point on the elliptical
+ # cross-section through the Earth towards the satellite, where the
+ # major axis is a and the minor is r. With the ellipse centered on the
+ # Earth and the satellite on the y-axis (at y = a + h = sat_dist), the
+ # equation for an ellipse and some calculus gives us the tangent point
+ # (x0, y0) as:
+ # y0 = a**2 / sat_dist
+ # x0 = r * np.sqrt(1 - a**2 / sat_dist**2)
+ # which gives:
+ # sin_c = x0 / np.hypot(x0, sat_dist - y0)
+ # tan_c = x0 / (sat_dist - y0)
+ # A bit of algebra combines these to give directly:
+ sin_c = r / np.sqrt(sat_dist ** 2 - a ** 2 + r ** 2)
+ tan_c = r / np.sqrt(sat_dist ** 2 - a ** 2)
+
+ # Using Napier's rules for right spherical triangles R2 and R6,
+ # (See https://en.wikipedia.org/wiki/Spherical_trigonometry), we can
+ # solve for arc angles b and a, which are our x and y scanning angles,
+ # respectively.
+ coords = np.vstack([np.arctan(np.cos(angleA) * tan_c), # R6
+ np.arcsin(np.sin(angleA) * sin_c)]) # R2
+
+ # Need to multiply scanning angles by satellite height to get to the
+ # actual native coordinates for the projection.
+ coords *= h
+ coords += np.array([[false_easting], [false_northing]])
+ return sgeom.LinearRing(coords.T)
+
+ @property
+ def boundary(self):
+ """Return boundary."""
+ return self._boundary
=====================================
pyresample/version.py
=====================================
@@ -26,9 +26,9 @@ def get_keywords() -> Dict[str, str]:
# 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 = " (tag: v1.28.4)"
- git_full = "8944215aab82ba2e941d4dfa4141ef23ac63ee2b"
- git_date = "2024-07-01 15:42:30 -0500"
+ git_refnames = " (HEAD -> main, tag: v1.29.0)"
+ git_full = "63b90b7964ab28001a370463de4f11c80e93a63e"
+ git_date = "2024-07-31 21:14:48 -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/-/commit/95525b45c770b679eeb20b01435795c0e1581785
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyresample/-/commit/95525b45c770b679eeb20b01435795c0e1581785
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/20240804/f32b39a5/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list