[Git][debian-gis-team/pyresample][upstream] New upstream version 1.12.3
Antonio Valentino
gitlab at salsa.debian.org
Wed Jul 10 07:23:53 BST 2019
Antonio Valentino pushed to branch upstream at Debian GIS Project / pyresample
Commits:
b3d4a977 by Antonio Valentino at 2019-07-10T05:56:19Z
New upstream version 1.12.3
- - - - -
16 changed files:
- CHANGELOG.md
- docs/source/geo_def.rst
- docs/source/geometry_utils.rst
- pyresample/_cartopy.py
- pyresample/_spatial_mp.py
- pyresample/area_config.py
- pyresample/geometry.py
- pyresample/kd_tree.py
- pyresample/test/__init__.py
- pyresample/test/test_files/areas.yaml
- pyresample/test/test_geometry.py
- + pyresample/test/test_spatial_mp.py
- pyresample/test/test_utils.py
- pyresample/utils/__init__.py
- pyresample/utils/_proj4.py
- pyresample/version.py
Changes:
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,51 @@
+## Version 1.12.3 (2019/05/17)
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 193](https://github.com/pytroll/pyresample/pull/193) - Fix striding slicing in AreaDefinition
+
+In this release 1 pull request was closed.
+
+
+## Version 1.12.2 (2019/05/10)
+
+### Issues Closed
+
+* [Issue 187](https://github.com/pytroll/pyresample/issues/187) - Numerous `RuntimeWarning`s when resampling
+
+In this release 1 issue was closed.
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 190](https://github.com/pytroll/pyresample/pull/190) - Fix aggregate method using non-serializable internal function
+* [PR 189](https://github.com/pytroll/pyresample/pull/189) - Fix dask race condition in KDTree resampling
+
+#### Features added
+
+* [PR 183](https://github.com/pytroll/pyresample/pull/183) - Fix bb computation to generate areas with equal h and v resolutions
+
+In this release 3 pull requests were closed.
+
+
+## Version 1.12.1 (2019/04/24)
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 186](https://github.com/pytroll/pyresample/pull/186) - Fix support for pyproj-2 EPSG syntax
+
+#### Documentation changes
+
+* [PR 185](https://github.com/pytroll/pyresample/pull/185) - Fix argument order in get_area_def doc
+
+In this release 2 pull requests were closed.
+
+
## Version 1.12.0 (2019/04/06)
### Issues Closed
@@ -15,8 +63,9 @@ In this release 1 issue was closed.
#### Features added
* [PR 182](https://github.com/pytroll/pyresample/pull/182) - Implement striding and aggregation for Swath- and AreaDefinitions
+* [PR 180](https://github.com/pytroll/pyresample/pull/180) - Remove radians from create_area_def and allow compatibility with pyproj-2.0+
-In this release 2 pull requests were closed.
+In this release 3 pull requests were closed.
## Version 1.11.2 (2019/03/18)
=====================================
docs/source/geo_def.rst
=====================================
@@ -51,11 +51,13 @@ where
* **upper_right_x**: projection x coordinate of upper right corner of upper right pixel
* **upper_right_y**: projection y coordinate of upper right corner of upper right pixel
-Below is an example of creating an ``AreaDefinition``:
+Below are three examples of creating an ``AreaDefinition``:
.. doctest::
>>> from pyresample.geometry import AreaDefinition
+
+ >>> # a) Using a projection dictionary
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
@@ -74,6 +76,37 @@ Below is an example of creating an ``AreaDefinition``:
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+ >>> # b) Using an explicit proj4 string
+ >>> proj_string = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
+ >>> area_def = AreaDefinition(area_id, description, proj_id, proj_string,
+ ... width, height, area_extent)
+ >>> print(area_def)
+ Area ID: ease_sh
+ Description: Antarctic EASE grid
+ Projection ID: ease_sh
+ Projection: {'a': '6371228.0', 'lat_0': '-90.0', 'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}
+ Number of columns: 425
+ Number of rows: 425
+ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+
+ >>> # c) Using an EPSG code in a proj4 string
+ >>> proj_string = '+init=EPSG:3409' # Use 'EPSG:3409' with pyproj 2.0+
+ >>> area_def = AreaDefinition(area_id, description, proj_id, proj_string,
+ ... width, height, area_extent)
+ >>> print(area_def)
+ Area ID: ease_sh
+ Description: Antarctic EASE grid
+ Projection ID: ease_sh
+ Projection: {'init': 'EPSG:3409'}
+ Number of columns: 425
+ Number of rows: 425
+ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
+
+.. note::
+
+ When using pyproj 2.0+, please use the new ``'EPSG:XXXX'`` syntax
+ as the old ``'+init=EPSG:XXXX'`` is no longer supported.
+
Creating an ``AreaDefinition`` can be complex if you don't know everything
about the region being described. Pyresample provides multiple utilities
for creating areas as well as storing them on disk for repeated use. See
=====================================
docs/source/geometry_utils.rst
=====================================
@@ -212,10 +212,10 @@ from_ul_corner
Loading from disk
-----------------
-The :func:`~pyresample.utils.load_area` function can be used to
+The :func:`~pyresample.area_config.load_area` function can be used to
parse area definitions from a configuration file by giving it the
-area file name and regions you wish to load. :func:`~pyresample.utils.load_area`
-takes advantage of :func:`~pyresample.utils.create_area_def`
+area file name and regions you wish to load. :func:`~pyresample.area_config.load_area`
+takes advantage of :func:`~pyresample.area_config.create_area_def`
and hence allows for the same arguments in the on-disk file.
Pyresample uses the YAML file format to store on-disk area definitions.
Below is an example YAML configuration file showing the various ways
@@ -344,6 +344,14 @@ an area might be specified.
resolution: 0.22542974631297721
units: deg
+ epsg:
+ area_id: ease_sh
+ description: Example of making an area definition using EPSG codes
+ projection:
+ init: EPSG:3410
+ shape: [425, 425]
+ area_extent: [-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625]
+
.. note::
The `lower_left_xy` and `upper_right_xy` items give the coordinates of the
@@ -351,13 +359,18 @@ an area might be specified.
projection coordinates are longitudes and latitudes, it is expected to
provide the extent in `longitude, latitude` order.
+.. note::
+
+ When using pyproj 2.0+, please use the new ``'EPSG: XXXX'`` syntax
+ as the old ``'init: EPSG:XXXX'`` is no longer supported.
+
If we assume the YAML content is stored in an ``areas.yaml`` file, we can
read a single ``AreaDefinition`` named ``corner`` by doing:
.. doctest::
- >>> from pyresample import utils
- >>> area_def = utils.load_area('areas.yaml', 'corner')
+ >>> from pyresample import load_area
+ >>> area_def = load_area('areas.yaml', 'corner')
>>> print(area_def)
Area ID: corner
Description: Example of making an area definition using shape, upper_left_extent, and resolution
@@ -371,7 +384,7 @@ series of arguments:
.. doctest::
- >>> corner, boundary = utils.load_area('areas.yaml', 'corner', 'boundary')
+ >>> corner, boundary = load_area('areas.yaml', 'corner', 'boundary')
>>> print(boundary)
Area ID: ease_sh
Description: Example of making an area definition using shape and area_extent
@@ -410,8 +423,8 @@ An area definition dict can be read using
.. doctest::
- >>> from pyresample import utils
- >>> area = utils.load_area('areas.cfg', 'ease_nh')
+ >>> from pyresample import load_area
+ >>> area = load_area('areas.cfg', 'ease_nh')
>>> print(area)
Area ID: ease_nh
Description: Arctic EASE grid
@@ -427,7 +440,7 @@ Several area definitions can be read at once using the region names in an argume
.. doctest::
- >>> nh_def, sh_def = utils.load_area('areas.cfg', 'ease_nh', 'ease_sh')
+ >>> nh_def, sh_def = load_area('areas.cfg', 'ease_nh', 'ease_sh')
>>> print(sh_def)
Area ID: ease_sh
Description: Antarctic EASE grid
=====================================
pyresample/_cartopy.py
=====================================
@@ -19,6 +19,8 @@
from logging import getLogger
import numpy as np
+import pyproj
+import warnings
try:
from xarray import DataArray
@@ -59,6 +61,9 @@ def _globe_from_proj4(proj4_terms):
class _PROJ4Projection(ccrs.Projection):
def __init__(self, proj4_terms, globe=None, bounds=None):
+ if 'EPSG' in proj4_terms.upper():
+ warnings.warn('Converting EPSG projection to proj4 string, which is a potentially lossy transformation')
+ proj4_terms = pyproj.Proj(proj4_terms).definition_string().strip()
terms = proj4_str_to_dict(proj4_terms)
globe = _globe_from_proj4(terms) if globe is None else globe
=====================================
pyresample/_spatial_mp.py
=====================================
@@ -24,7 +24,7 @@ import pyproj
import multiprocessing as mp
import warnings
-from pyresample.utils import proj4_str_to_dict
+from pyresample.utils import proj4_str_to_dict, is_pyproj2
try:
@@ -112,8 +112,12 @@ class cKDTree_MP(object):
class BaseProj(pyproj.Proj):
def __init__(self, projparams=None, preserve_units=True, **kwargs):
+ # Copy dict-type arguments as they will be modified in-place
+ if isinstance(projparams, dict):
+ projparams = projparams.copy()
+
# Pyproj<2 uses __new__ to initiate data and does not define its own __init__ method.
- if pyproj.__version__ >= '2':
+ if is_pyproj2():
# If init is found in any of the data, override any other area parameters.
if 'init' in kwargs:
warnings.warn('init="EPSG:XXXX" is no longer supported. Use "EPSG:XXXX" as a proj string instead')
@@ -127,10 +131,18 @@ class BaseProj(pyproj.Proj):
projparams = proj4_str_to_dict(projparams)
warnings.warn(warn_msg)
projparams = projparams.pop('init')
+ # New syntax 'EPSG:XXXX'
+ if 'EPSG' in kwargs or (isinstance(projparams, dict) and 'EPSG' in projparams):
+ if 'EPSG' in kwargs:
+ epsg_code = kwargs.pop('EPSG')
+ else:
+ epsg_code = projparams.pop('EPSG')
+ projparams = 'EPSG:{}'.format(epsg_code)
+
super(BaseProj, self).__init__(projparams=projparams, preserve_units=preserve_units, **kwargs)
def is_latlong(self):
- if pyproj.__version__ >= '2':
+ if is_pyproj2():
return self.crs.is_geographic
return super(BaseProj, self).is_latlong()
=====================================
pyresample/area_config.py
=====================================
@@ -305,10 +305,10 @@ def get_area_def(area_id, area_name, proj_id, proj4_args, width, height, area_ex
-----------
area_id : str
ID of area
- proj_id : str
- ID of projection
area_name :str
Description of area
+ proj_id : str
+ ID of projection
proj4_args : list, dict, or str
Proj4 arguments as list of arguments or string
width : int
@@ -453,9 +453,12 @@ def create_area_def(area_id, projection, width=None, height=None, area_extent=No
# Fills in missing information to attempt to create an area definition.
if area_extent is None or shape is None:
- area_extent, shape = _extrapolate_information(area_extent, shape, center, radius, resolution,
- upper_left_extent, units, p, proj_dict)
- return _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **kwargs)
+ area_extent, shape, resolution = \
+ _extrapolate_information(area_extent, shape, center, radius,
+ resolution, upper_left_extent, units,
+ p, proj_dict)
+ return _make_area(area_id, description, proj_id, proj_dict, shape,
+ area_extent, resolution=resolution, **kwargs)
def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **kwargs):
@@ -464,6 +467,7 @@ def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **k
from pyresample.geometry import DynamicAreaDefinition
# Remove arguments that are only for DynamicAreaDefinition.
optimize_projection = kwargs.pop('optimize_projection', False)
+ resolution = kwargs.pop('resolution', None)
# If enough data is provided, create an AreaDefinition. If only shape or area_extent are found, make a
# DynamicAreaDefinition. If not enough information was provided, raise a ValueError.
height, width = (None, None)
@@ -474,7 +478,7 @@ def _make_area(area_id, description, proj_id, proj_dict, shape, area_extent, **k
return DynamicAreaDefinition(area_id=area_id, description=description, projection=proj_dict, width=width,
height=height, area_extent=area_extent, rotation=kwargs.get('rotation'),
- optimize_projection=optimize_projection)
+ resolution=resolution, optimize_projection=optimize_projection)
def _get_proj_data(projection):
@@ -523,7 +527,7 @@ def _distance_from_center_forward(var, center, p):
# Since the distance between longitudes and latitudes is not constant in
# most projections, there must be reference point to start from.
if center is None:
- raise ValueError('center must be given to convert radius or resolution from an angle to meters')
+ center = (0, 0)
center_as_angle = p(*center, inverse=True, errcheck=True)
pole = 90
@@ -700,7 +704,7 @@ def _extrapolate_information(area_extent, shape, center, radius, resolution, upp
upper_left_extent[0], upper_left_extent[1] - 2 * radius[1], upper_left_extent[0] + 2 * radius[0],
upper_left_extent[1])
area_extent = _validate_variable(area_extent, new_area_extent, 'area_extent', ['upper_left_extent', 'radius'])
- return area_extent, shape
+ return area_extent, shape, resolution
def _format_list(var, name):
=====================================
pyresample/geometry.py
=====================================
@@ -31,7 +31,7 @@ from logging import getLogger
import numpy as np
import yaml
-from pyproj import Geod
+from pyproj import Geod, transform
from pyresample import CHUNK_SIZE, utils
from pyresample._spatial_mp import Cartesian, Cartesian_MP, Proj, Proj_MP
@@ -551,6 +551,12 @@ class SwathDefinition(CoordinateDefinition):
"""Copy the current swath."""
return SwathDefinition(self.lons, self.lats)
+ @staticmethod
+ def _do_transform(src, dst, lons, lats, alt):
+ """Helper for 'aggregate' method."""
+ x, y, z = transform(src, dst, lons, lats, alt)
+ return np.dstack((x, y, z))
+
def aggregate(self, **dims):
"""Aggregate the current swath definition by averaging.
@@ -560,20 +566,15 @@ class SwathDefinition(CoordinateDefinition):
import pyproj
import dask.array as da
- def do_transform(src, dst, lons, lats, alt):
- import pyproj
- x, y, z = pyproj.transform(src, dst, lons, lats, alt)
- return np.dstack((x, y, z))
-
geocent = pyproj.Proj(proj='geocent')
latlong = pyproj.Proj(proj='latlong')
- res = da.map_blocks(do_transform, latlong, geocent,
+ res = da.map_blocks(self._do_transform, latlong, geocent,
self.lons.data, self.lats.data,
da.zeros_like(self.lons.data), new_axis=[2],
chunks=(self.lons.chunks[0], self.lons.chunks[1], 3))
res = DataArray(res, dims=['y', 'x', 'coord'], coords=self.lons.coords)
res = res.coarsen(**dims).mean()
- lonlatalt = da.map_blocks(do_transform, geocent, latlong,
+ lonlatalt = da.map_blocks(self._do_transform, geocent, latlong,
res[:, :, 0].data, res[:, :, 1].data,
res[:, :, 2].data, new_axis=[2],
chunks=res.data.chunks)
@@ -676,6 +677,33 @@ class SwathDefinition(CoordinateDefinition):
new_proj.update(proj_dict)
return new_proj
+ def _compute_uniform_shape(self):
+ """Compute the height and width of a domain to have uniform resolution across dimensions."""
+ g = Geod(ellps='WGS84')
+ leftlons = self.lons[:, 0]
+ leftlons = leftlons.where(leftlons.notnull(), drop=True)
+ rightlons = self.lons[:, -1]
+ rightlons = rightlons.where(rightlons.notnull(), drop=True)
+ middlelons = self.lons[:, int(self.lons.shape[1] / 2)]
+ middlelons = middlelons.where(middlelons.notnull(), drop=True)
+ leftlats = self.lats[:, 0]
+ leftlats = leftlats.where(leftlats.notnull(), drop=True)
+ rightlats = self.lats[:, -1]
+ rightlats = rightlats.where(rightlats.notnull(), drop=True)
+ middlelats = self.lats[:, int(self.lats.shape[1] / 2)]
+ middlelats = middlelats.where(middlelats.notnull(), drop=True)
+
+ az1, az2, width1 = g.inv(leftlons[0], leftlats[0], rightlons[0], rightlats[0])
+ 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)
+ width = int(width * 1.1 / resolution)
+ height = int(height * 1.1 / resolution)
+ return height, width
+
def compute_optimal_bb_area(self, proj_dict=None):
"""Compute the "best" bounding box area for this swath with `proj_dict`.
@@ -683,16 +711,16 @@ class SwathDefinition(CoordinateDefinition):
which case the right projection angle `alpha` is computed from the
swath centerline. For other projections, only the appropriate center of
projection and area extents are computed.
+
+ The height and width are computed so that the resolution is
+ approximately the same across dimensions.
"""
if proj_dict is None:
proj_dict = {}
projection = proj_dict.setdefault('proj', 'omerc')
area_id = projection + '_otf'
description = 'On-the-fly ' + projection + ' area'
- lines, cols = self.lons.shape
- width = int(cols * 1.1)
- height = int(lines * 1.1)
-
+ height, width = self._compute_uniform_shape()
proj_dict = self.compute_bb_proj_params(proj_dict)
area = DynamicAreaDefinition(area_id, description, proj_dict)
@@ -768,10 +796,10 @@ class DynamicAreaDefinition(object):
x_resolution, y_resolution = resolution
except TypeError:
x_resolution = y_resolution = resolution
- width = int(np.rint((corners[2] - corners[0]) * 1.0 /
- x_resolution + 1))
- height = int(np.rint((corners[3] - corners[1]) * 1.0 /
- y_resolution + 1))
+ width = int(np.rint((corners[2] - corners[0]) * 1.0
+ / x_resolution + 1))
+ height = int(np.rint((corners[3] - corners[1]) * 1.0
+ / y_resolution + 1))
area_extent = (corners[0] - x_resolution / 2,
corners[1] - y_resolution / 2,
@@ -1712,8 +1740,10 @@ class AreaDefinition(BaseDefinition):
yslice, xslice = key
# Get actual values, replace Nones
yindices = yslice.indices(self.height)
+ total_rows = int((yindices[1] - yindices[0]) / yindices[2])
ystopactual = yindices[1] - (yindices[1] - 1) % yindices[2]
xindices = xslice.indices(self.width)
+ total_cols = int((xindices[1] - xindices[0]) / xindices[2])
xstopactual = xindices[1] - (xindices[1] - 1) % xindices[2]
yslice = slice(yindices[0], ystopactual, yindices[2])
xslice = slice(xindices[0], xstopactual, xindices[2])
@@ -1725,8 +1755,8 @@ class AreaDefinition(BaseDefinition):
new_area = AreaDefinition(self.area_id, self.description,
self.proj_id, self.proj_dict,
- int((xslice.stop - xslice.start) / xslice.step),
- int((yslice.stop - yslice.start) / yslice.step),
+ total_cols,
+ total_rows,
new_area_extent)
new_area.crop_offset = (self.crop_offset[0] + yslice.start,
self.crop_offset[1] + xslice.start)
=====================================
pyresample/kd_tree.py
=====================================
@@ -892,13 +892,13 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
def lonlat2xyz(lons, lats):
-
R = 6370997.0
- x_coords = R * da.cos(da.deg2rad(lats)) * da.cos(da.deg2rad(lons))
- y_coords = R * da.cos(da.deg2rad(lats)) * da.sin(da.deg2rad(lons))
- z_coords = R * da.sin(da.deg2rad(lats))
+ x_coords = R * np.cos(np.deg2rad(lats)) * np.cos(np.deg2rad(lons))
+ y_coords = R * np.cos(np.deg2rad(lats)) * np.sin(np.deg2rad(lons))
+ z_coords = R * np.sin(np.deg2rad(lats))
- return da.stack(
+ stack = np.stack if isinstance(lons, np.ndarray) else da.stack
+ return stack(
(x_coords.ravel(), y_coords.ravel(), z_coords.ravel()), axis=-1)
@@ -922,7 +922,7 @@ def query_no_distance(target_lons, target_lats, valid_output_index,
coords = lonlat2xyz(target_lons_valid, target_lats_valid)
distance_array, index_array = kdtree.query(
- coords.compute(),
+ coords,
k=neighbours,
eps=epsilon,
distance_upper_bound=radius,
@@ -946,6 +946,20 @@ def query_no_distance(target_lons, target_lats, valid_output_index,
return res_ia
+def _my_index(index_arr, vii, data_arr, vii_slices=None, ia_slices=None,
+ fill_value=np.nan):
+ """Helper function for 'get_sample_from_neighbour_info'."""
+ vii_slices = tuple(
+ x if x is not None else vii.ravel() for x in vii_slices)
+ mask_slices = tuple(
+ x if x is not None else (index_arr == -1) for x in ia_slices)
+ ia_slices = tuple(
+ x if x is not None else index_arr for x in ia_slices)
+ res = data_arr[vii_slices][ia_slices]
+ res[mask_slices] = fill_value
+ return res
+
+
class XArrayResamplerNN(object):
def __init__(self,
source_geo_def,
@@ -1176,18 +1190,6 @@ class XArrayResamplerNN(object):
# FUTURE: when we allow more than one neighbor
# neighbors_dim = i + 3
- def _my_index(index_arr, vii, data_arr, vii_slices=None,
- ia_slices=None, fill_value=np.nan):
- vii_slices = tuple(
- x if x is not None else vii.ravel() for x in vii_slices)
- mask_slices = tuple(
- x if x is not None else (index_arr == -1) for x in ia_slices)
- ia_slices = tuple(
- x if x is not None else index_arr for x in ia_slices)
- res = data_arr[vii_slices][ia_slices]
- res[mask_slices] = fill_value
- return res
-
new_data = data.data.reshape(flat_src_shape)
vii = vii.ravel()
dst_adims = [dst_dim_to_ind[dim] for dim in dst_dims]
=====================================
pyresample/test/__init__.py
=====================================
@@ -37,6 +37,7 @@ from pyresample.test import (
test_ewa_fornav,
test_bilinear,
test_data_reduce,
+ test_spatial_mp
)
import unittest
@@ -59,6 +60,7 @@ def suite():
mysuite.addTests(test_ewa_fornav.suite())
mysuite.addTests(test_bilinear.suite())
mysuite.addTests(test_data_reduce.suite())
+ mysuite.addTests(test_spatial_mp.suite())
return mysuite
=====================================
pyresample/test/test_files/areas.yaml
=====================================
@@ -218,4 +218,26 @@ test_latlong:
width: 3473
area_extent:
lower_left_xy: [-0.08115781021773638, 0.4038691889114878]
- upper_right_xy: [0.08115781021773638, 0.5427973973702365]
\ No newline at end of file
+ upper_right_xy: [0.08115781021773638, 0.5427973973702365]
+
+test_dynamic_resolution:
+ description: Dynamic with resolution specified in meters
+ projection:
+ proj: lcc
+ lon_0: -95.0
+ lat_0: 25.0
+ lat_1: 25.0
+ ellps: WGS84
+ resolution: [1000, 1000]
+
+test_dynamic_resolution_ll:
+ description: Dynamic with resolution specified in degrees
+ projection:
+ proj: longlat
+ lat_0: 27.12
+ lon_0: -81.36
+ ellps: WGS84
+ resolution:
+ dy: 1
+ dx: 1
+ units: deg
=====================================
pyresample/test/test_geometry.py
=====================================
@@ -6,7 +6,7 @@ import sys
import numpy as np
-from pyresample import geo_filter, geometry
+from pyresample import geo_filter, geometry, parse_area_file
from pyresample.geometry import (IncompatibleAreas,
combine_area_extents_vertical,
concatenate_area_defs)
@@ -70,6 +70,8 @@ class Test(unittest.TestCase):
def test_cartopy_crs(self):
"""Test conversion from area definition to cartopy crs"""
+ from pyresample import utils
+
europe = geometry.AreaDefinition(area_id='areaD',
description='Europe (3km, HRV, VTC)',
proj_id='areaD',
@@ -111,7 +113,28 @@ class Test(unittest.TestCase):
np.fabs(area_def.area_extent[3] - area_def.area_extent[1])) / 100.
self.assertEqual(crs.threshold, thresh_exp)
+ # EPSG projection
+ projections = ['+init=EPSG:6932']
+ if utils.is_pyproj2():
+ projections.append('EPSG:6932')
+
+ for projection in projections:
+ area = geometry.AreaDefinition(
+ area_id='ease-sh-2.0',
+ description='25km EASE Grid 2.0 (Southern Hemisphere)',
+ proj_id='ease-sh-2.0',
+ projection=projection,
+ width=123, height=123,
+ area_extent=[-40000., -40000., 40000., 40000.])
+ with patch('pyresample._cartopy.warnings.warn') as warn:
+ # Test that user warning has been issued (EPSG to proj4 string is potentially lossy)
+ area.to_cartopy_crs()
+ warn.assert_called()
+
def test_create_areas_def(self):
+ from pyresample import utils
+ import yaml
+
area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
'areaD',
{'a': '6378144.0',
@@ -126,18 +149,93 @@ class Test(unittest.TestCase):
-909968.64000000001,
1029087.28,
1490031.3600000001])
- import yaml
- res = yaml.load(area_def.create_areas_def())
- expected = yaml.load(('areaD:\n description: Europe (3km, HRV, VTC)\n'
- ' projection:\n a: 6378144.0\n b: 6356759.0\n'
- ' lat_0: 50.0\n lat_ts: 50.0\n lon_0: 8.0\n'
- ' proj: stere\n shape:\n height: 800\n'
- ' width: 800\n area_extent:\n'
- ' lower_left_xy: [-1370912.72, -909968.64]\n'
- ' upper_right_xy: [1029087.28, 1490031.36]\n'))
+ res = yaml.safe_load(area_def.create_areas_def())
+ expected = yaml.safe_load(('areaD:\n description: Europe (3km, HRV, VTC)\n'
+ ' projection:\n a: 6378144.0\n b: 6356759.0\n'
+ ' lat_0: 50.0\n lat_ts: 50.0\n lon_0: 8.0\n'
+ ' proj: stere\n shape:\n height: 800\n'
+ ' width: 800\n area_extent:\n'
+ ' lower_left_xy: [-1370912.72, -909968.64]\n'
+ ' upper_right_xy: [1029087.28, 1490031.36]\n'))
self.assertDictEqual(res, expected)
+ # EPSG
+ projections = {'+init=epsg:3006': 'init: epsg:3006'}
+ if utils.is_pyproj2():
+ projections['EPSG:3006'] = 'EPSG: 3006'
+
+ for projection, epsg_yaml in projections.items():
+ area_def = geometry.AreaDefinition('baws300_sweref99tm', 'BAWS, 300m resolution, sweref99tm',
+ 'sweref99tm',
+ projection,
+ 4667,
+ 4667,
+ [-49739, 5954123, 1350361, 7354223])
+ res = yaml.safe_load(area_def.create_areas_def())
+ expected = yaml.safe_load(('baws300_sweref99tm:\n'
+ ' description: BAWS, 300m resolution, sweref99tm\n'
+ ' projection:\n'
+ ' {epsg}\n'
+ ' shape:\n'
+ ' height: 4667\n'
+ ' width: 4667\n'
+ ' area_extent:\n'
+ ' lower_left_xy: [-49739, 5954123]\n'
+ ' upper_right_xy: [1350361, 7354223]'.format(epsg=epsg_yaml)))
+ self.assertDictEqual(res, expected)
+
+ def test_parse_area_file(self):
+ from pyresample import utils
+
+ expected = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
+ 'areaD',
+ {'a': '6378144.0',
+ 'b': '6356759.0',
+ 'lat_0': '50.00',
+ 'lat_ts': '50.00',
+ 'lon_0': '8.00',
+ 'proj': 'stere'},
+ 800,
+ 800,
+ [-1370912.72,
+ -909968.64000000001,
+ 1029087.28,
+ 1490031.3600000001])
+ yaml_str = ('areaD:\n description: Europe (3km, HRV, VTC)\n'
+ ' projection:\n a: 6378144.0\n b: 6356759.0\n'
+ ' lat_0: 50.0\n lat_ts: 50.0\n lon_0: 8.0\n'
+ ' proj: stere\n shape:\n height: 800\n'
+ ' width: 800\n area_extent:\n'
+ ' lower_left_xy: [-1370912.72, -909968.64]\n'
+ ' upper_right_xy: [1029087.28, 1490031.36]\n')
+ area_def = parse_area_file(yaml_str, 'areaD')[0]
+ self.assertEqual(area_def, expected)
+
+ # EPSG
+ projections = {'+init=epsg:3006': 'init: epsg:3006'}
+ if utils.is_pyproj2():
+ projections['EPSG:3006'] = 'EPSG: 3006'
+ for projection, epsg_yaml in projections.items():
+ expected = geometry.AreaDefinition('baws300_sweref99tm', 'BAWS, 300m resolution, sweref99tm',
+ 'sweref99tm',
+ projection,
+ 4667,
+ 4667,
+ [-49739, 5954123, 1350361, 7354223])
+ yaml_str = ('baws300_sweref99tm:\n'
+ ' description: BAWS, 300m resolution, sweref99tm\n'
+ ' projection:\n'
+ ' {epsg}\n'
+ ' shape:\n'
+ ' height: 4667\n'
+ ' width: 4667\n'
+ ' area_extent:\n'
+ ' lower_left_xy: [-49739, 5954123]\n'
+ ' upper_right_xy: [1350361, 7354223]'.format(epsg=epsg_yaml))
+ area_def = parse_area_file(yaml_str, 'baws300_sweref99tm')[0]
+ self.assertEqual(area_def, expected)
+
def test_base_type(self):
lons1 = np.arange(-135., +135, 50.)
lats = np.ones_like(lons1) * 70.
@@ -916,16 +1014,21 @@ class Test(unittest.TestCase):
self.assertEqual(slice_y, slice(158, 515, None))
# totally different area
- area_to_cover = geometry.AreaDefinition('epsg4326', 'Global equal latitude/longitude grid for global sphere',
- 'epsg4326',
- {"init": 'EPSG:4326'},
- 8192,
- 4096,
- [-180.0, -90.0, 180.0, 90.0])
-
- slice_x, slice_y = area_def.get_area_slices(area_to_cover)
- self.assertEqual(slice_x, slice(46, 3667, None))
- self.assertEqual(slice_y, slice(52, 3663, None))
+ projections = [{"init": 'EPSG:4326'}]
+ if utils.is_pyproj2():
+ projections.append('EPSG:4326')
+ for projection in projections:
+ area_to_cover = geometry.AreaDefinition(
+ 'epsg4326', 'Global equal latitude/longitude grid for global sphere',
+ 'epsg4326',
+ projection,
+ 8192,
+ 4096,
+ [-180.0, -90.0, 180.0, 90.0])
+
+ slice_x, slice_y = area_def.get_area_slices(area_to_cover)
+ self.assertEqual(slice_x, slice(46, 3667, None))
+ self.assertEqual(slice_y, slice(52, 3663, None))
def test_get_area_slices_nongeos(self):
"""Check area slicing for non-geos projections."""
@@ -964,6 +1067,8 @@ class Test(unittest.TestCase):
def test_proj_str(self):
from collections import OrderedDict
+ from pyresample import utils
+
proj_dict = OrderedDict()
proj_dict['proj'] = 'stere'
proj_dict['a'] = 6378144.0
@@ -985,6 +1090,20 @@ class Test(unittest.TestCase):
self.assertEqual(area.proj_str,
'+a=6378144.0 +b=6356759.0 +lat_0=50.0 +lat_ts=50.0 +lon_0=8.0 +no_rot +proj=stere')
+ # EPSG
+ projections = ['+init=EPSG:6932']
+ if utils.is_pyproj2():
+ projections.append('EPSG:6932')
+ for projection in projections:
+ area = geometry.AreaDefinition(
+ area_id='ease-sh-2.0',
+ description='25km EASE Grid 2.0 (Southern Hemisphere)',
+ proj_id='ease-sh-2.0',
+ projection=projection,
+ width=123, height=123,
+ area_extent=[-40000., -40000., 40000., 40000.])
+ self.assertEqual(area.proj_str, projection)
+
def test_striding(self):
"""Test striding AreaDefinitions."""
from pyresample import utils
@@ -1009,6 +1128,7 @@ class Test(unittest.TestCase):
area_extent[1] + 3 * area_def.pixel_size_y,
area_extent[2] - 3 * area_def.pixel_size_x,
area_extent[3]))
+ self.assertEqual(reduced_area.shape, (928, 928))
def assert_np_dict_allclose(dict1, dict2):
@@ -1247,25 +1367,27 @@ class TestSwathDefinition(unittest.TestCase):
def test_compute_optimal_bb(self):
"""Test computing the bb area."""
+ import xarray as xr
lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
[80.84000396728516, 60.74200439453125, 34.08500289916992],
[67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
-
+ lats = xr.DataArray(lats)
lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
[79.11000061035156, 7.284000396728516, -5.107000350952148],
[81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+ lons = xr.DataArray(lons)
area = geometry.SwathDefinition(lons, lats)
res = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'})
- np.testing.assert_allclose(res.area_extent, [-2348379.728104, 2284625.526467,
- 2432121.058435, 11719235.223912])
+ np.testing.assert_allclose(res.area_extent, [-2348379.728104, 3228086.496211,
+ 2432121.058435, 10775774.254169])
proj_dict = {'gamma': 0.0, 'lonc': -11.391744043133668,
'ellps': 'WGS84', 'proj': 'omerc',
'alpha': 9.185764390923012, 'lat_0': -0.2821013754097188}
assert_np_dict_allclose(res.proj_dict, proj_dict)
- self.assertEqual(res.shape, (3, 3))
+ self.assertEqual(res.shape, (6, 3))
def test_aggregation(self):
"""Test aggregation on SwathDefinitions."""
@@ -1472,11 +1594,16 @@ class TestStackedAreaDefinition(unittest.TestCase):
from pyresample.geometry import DynamicAreaDefinition
from pyresample.area_config import DataArray
from pyresample.area_config import create_area_def as cad
+ from pyresample import utils
+ import pyproj
area_id = 'ease_sh'
description = 'Antarctic EASE grid'
projection_list = [{'proj': 'laea', 'lat_0': -90, 'lon_0': 0, 'a': 6371228.0, 'units': 'm'},
- '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m']
+ '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m',
+ '+init=EPSG:3409']
+ if utils.is_pyproj2():
+ projection_list.append('EPSG:3409')
proj_id = 'ease_sh'
shape = (425, 850)
upper_left_extent = (-5326849.0625, 5326849.0625)
@@ -1511,7 +1638,7 @@ class TestStackedAreaDefinition(unittest.TestCase):
radius=essentials[1], description=description, units=units, rotation=45))
except ValueError:
pass
- self.assertEqual(len(area_list), 4)
+ self.assertEqual(len(area_list), 8 if utils.is_pyproj2() else 6)
# Tests that specifying units through xarrays works.
area_list.append(cad(area_id, projection_list[1], shape=shape,
@@ -1544,6 +1671,18 @@ class TestStackedAreaDefinition(unittest.TestCase):
self.assertEqual(area_def.shape, (101, 90))
# Checks every area definition made
for area_def in area_list:
+ if 'EPSG' in area_def.proj_dict or 'init' in area_def.proj_dict:
+ # Use formal definition of EPSG projections to make them comparable to the base definition
+ proj_def = pyproj.Proj(area_def.proj_str).definition_string().strip()
+ area_def = area_def.copy(projection=proj_def)
+
+ # Remove extra attributes from the formal definition
+ if 'R' in area_def.proj_dict:
+ # pyproj < 2
+ area_def.proj_dict['a'] = area_def.proj_dict.pop('R')
+ for key in ['x_0', 'y_0', 'no_defs', 'b', 'init']:
+ area_def.proj_dict.pop(key, None)
+
self.assertEqual(area_def, base_def)
# Makes sure if shape or area_extent is found/given, a DynamicAreaDefinition is made.
@@ -1588,14 +1727,15 @@ class TestDynamicAreaDefinition(unittest.TestCase):
lats = [[66, 67, 68, 69.],
[58, 59, 60, 61],
[50, 51, 52, 53]]
- sdef = geometry.SwathDefinition(lons, lats)
+ import xarray as xr
+ sdef = geometry.SwathDefinition(xr.DataArray(lons), xr.DataArray(lats))
result = area.freeze(sdef,
resolution=1000)
np.testing.assert_allclose(result.area_extent,
- [-336277.698941, 5047207.008079,
- 192456.651909, 8215588.023806])
+ [-336277.698941, 5513145.392745,
+ 192456.651909, 7749649.63914])
self.assertEqual(result.x_size, 4)
- self.assertEqual(result.y_size, 3)
+ self.assertEqual(result.y_size, 18)
def test_compute_domain(self):
"""Test computing size and area extent."""
=====================================
pyresample/test/test_spatial_mp.py
=====================================
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# pyresample, Resampling of remote sensing image data in python
+#
+# Copyright (C) 2014-2019 PyTroll developers
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""Testing the _spatial_mp module"""
+
+
+try:
+ from unittest import mock
+except ImportError:
+ # separate mock package py<3.3
+ import mock
+import pyproj
+import unittest
+
+from pyresample._spatial_mp import BaseProj
+
+
+class SpatialMPTest(unittest.TestCase):
+ @mock.patch('pyresample._spatial_mp.pyproj.Proj.__init__', return_value=None)
+ def test_base_proj_epsg(self, proj_init):
+ """Test Proj creation with EPSG codes"""
+ if pyproj.__version__ < '2':
+ return self.skipTest(reason='pyproj 2+ only')
+
+ args = [
+ [None, {'init': 'EPSG:6932'}],
+ [{'init': 'EPSG:6932'}, {}],
+ [None, {'EPSG': '6932'}],
+ [{'EPSG': '6932'}, {}]
+ ]
+ for projparams, kwargs in args:
+ BaseProj(projparams, **kwargs)
+ proj_init.assert_called_with(projparams='EPSG:6932', preserve_units=mock.ANY)
+ proj_init.reset_mock()
+
+
+def suite():
+ """The test suite.
+ """
+ loader = unittest.TestLoader()
+ mysuite = unittest.TestSuite()
+ mysuite.addTest(loader.loadTestsFromTestCase(SpatialMPTest))
+ return mysuite
+
+
+if __name__ == '__main__':
+ unittest.main()
=====================================
pyresample/test/test_utils.py
=====================================
@@ -136,6 +136,27 @@ Number of rows: 4058
Area extent: (-0.0812, 0.4039, 0.0812, 0.5428)"""
self.assertEqual(test_latlong.__str__(), latlong_str)
+ def test_dynamic_area_parser_yaml(self):
+ """Test YAML area parser on dynamic areas."""
+ from pyresample import parse_area_file
+ from pyresample.geometry import DynamicAreaDefinition
+ test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area = parse_area_file(test_area_file, 'test_dynamic_resolution')[0]
+
+ self.assertIsInstance(test_area, DynamicAreaDefinition)
+ self.assertTrue(hasattr(test_area, 'resolution'))
+ self.assertEqual(test_area.resolution, (1000.0, 1000.0))
+
+ # lat/lon
+ from pyresample import parse_area_file
+ from pyresample.geometry import DynamicAreaDefinition
+ test_area_file = os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml')
+ test_area = parse_area_file(test_area_file, 'test_dynamic_resolution_ll')[0]
+
+ self.assertIsInstance(test_area, DynamicAreaDefinition)
+ self.assertTrue(hasattr(test_area, 'resolution'))
+ self.assertEqual(test_area.resolution, (1.0, 1.0))
+
def test_multiple_file_content(self):
from pyresample import parse_area_file
area_list = ["""ease_sh:
@@ -291,8 +312,23 @@ class TestMisc(unittest.TestCase):
np.testing.assert_almost_equal(a, 6378137.)
np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
+ def test_convert_proj_floats(self):
+ from collections import OrderedDict
+ import pyresample.utils as utils
+
+ pairs = [('proj', 'lcc'), ('ellps', 'WGS84'), ('lon_0', '-95'), ('no_defs', True)]
+ expected = OrderedDict([('proj', 'lcc'), ('ellps', 'WGS84'), ('lon_0', -95.0), ('no_defs', True)])
+ self.assertDictEqual(utils._proj4.convert_proj_floats(pairs), expected)
+
+ # EPSG
+ pairs = [('init', 'EPSG:4326'), ('EPSG', 4326)]
+ for pair in pairs:
+ expected = OrderedDict([pair])
+ self.assertDictEqual(utils._proj4.convert_proj_floats([pair]), expected)
+
def test_proj4_str_dict_conversion(self):
from pyresample import utils
+
proj_str = "+proj=lcc +ellps=WGS84 +lon_0=-95 +no_defs"
proj_dict = utils._proj4.proj4_str_to_dict(proj_str)
proj_str2 = utils._proj4.proj4_dict_to_str(proj_dict)
@@ -301,6 +337,18 @@ class TestMisc(unittest.TestCase):
self.assertIsInstance(proj_dict['lon_0'], float)
self.assertIsInstance(proj_dict2['lon_0'], float)
+ # EPSG
+ expected = {'+init=EPSG:4326': {'init': 'EPSG:4326'},
+ 'EPSG:4326': {'EPSG': 4326}}
+
+ for proj_str, proj_dict_exp in expected.items():
+ proj_dict = utils._proj4.proj4_str_to_dict(proj_str)
+ self.assertEqual(proj_dict, proj_dict_exp)
+ self.assertEqual(utils._proj4.proj4_dict_to_str(proj_dict), proj_str) # round-trip
+
+ # Invalid EPSG code (pyproj-2 syntax only)
+ self.assertRaises(ValueError, utils._proj4.proj4_str_to_dict, 'EPSG:XXXX')
+
def test_def2yaml_converter(self):
from pyresample import parse_area_file, convert_def_to_yaml
import tempfile
=====================================
pyresample/utils/__init__.py
=====================================
@@ -18,6 +18,7 @@
"""Miscellaneous utility functions for pyresample."""
from collections import Mapping
import numpy as np
+import pyproj
import warnings
from ._proj4 import (proj4_dict_to_str, proj4_str_to_dict, convert_proj_floats, proj4_radius_parameters) # noqa
@@ -234,3 +235,8 @@ def recursive_dict_update(d, u):
else:
d[k] = u[k]
return d
+
+
+def is_pyproj2():
+ """Determine whether the current pyproj version is >= 2.0"""
+ return pyproj.__version__ >= '2'
=====================================
pyresample/utils/_proj4.py
=====================================
@@ -16,6 +16,7 @@
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
from collections import OrderedDict
+import six
def convert_proj_floats(proj_pairs):
@@ -25,6 +26,9 @@ def convert_proj_floats(proj_pairs):
if len(x) == 1 or x[1] is True:
proj_dict[x[0]] = True
continue
+ if x[0] == 'EPSG':
+ proj_dict[x[0]] = x[1]
+ continue
try:
proj_dict[x[0]] = float(x[1])
@@ -39,6 +43,13 @@ def proj4_str_to_dict(proj4_str):
Note: Key only parameters will be assigned a value of `True`.
"""
+ if proj4_str.startswith('EPSG:'):
+ try:
+ code = int(proj4_str.split(':', 1)[1])
+ except ValueError as err:
+ six.raise_from(ValueError("Invalid EPSG code '{}': {}".format(proj4_str, err)),
+ None) # Suppresses original exception context in python 3
+ return OrderedDict(EPSG=code)
pairs = (x.split('=', 1) for x in proj4_str.replace('+', '').split(" "))
return convert_proj_floats(pairs)
@@ -50,6 +61,11 @@ def proj4_dict_to_str(proj4_dict, sort=False):
items = sorted(items)
params = []
for key, val in items:
+ if key == 'EPSG':
+ # If EPSG code is present, ignore other parameters
+ params = ['EPSG:{}'.format(val)]
+ break
+
key = str(key) if key.startswith('+') else '+' + str(key)
if key in ['+no_defs', '+no_off', '+no_rot']:
param = key
=====================================
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 = " (tag: v1.12.0)"
- git_full = "5c4480e87898aaed2950e775324fec0bd408853c"
- git_date = "2019-04-06 20:47:48 -0500"
+ git_refnames = " (tag: v1.12.3)"
+ git_full = "6b4367ce09ca6040c6107d2fafc31524eb45d1f8"
+ git_date = "2019-05-17 10:58:40 -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/b3d4a977d54561b8243577637786eb8ebda0500d
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pyresample/commit/b3d4a977d54561b8243577637786eb8ebda0500d
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/20190710/77c1d7be/attachment-0001.html>
More information about the Pkg-grass-devel
mailing list