[pyresample] 01/08: New upstream version 1.6.1
Antonio Valentino
a_valentino-guest at moszumanska.debian.org
Sat Sep 23 17:15:52 UTC 2017
This is an automated email from the git hooks/post-receive script.
a_valentino-guest pushed a commit to branch master
in repository pyresample.
commit 34c33ee8172d2f1084f3eafef3fac431d9b8bc79
Author: Antonio Valentino <antonio.valentino at tiscali.it>
Date: Sat Sep 23 17:06:01 2017 +0000
New upstream version 1.6.1
---
.bumpversion.cfg | 2 +-
changelog.rst | 50 +++++
pyresample/geometry.py | 171 +++++++++++++++++
pyresample/kd_tree.py | 80 ++++----
pyresample/test/test_ewa_ll2cr.py | 23 +--
pyresample/test/test_geometry.py | 385 ++++++++++++++++++++++++++++----------
pyresample/test/test_utils.py | 87 ++++++++-
pyresample/test/utils.py | 22 +++
pyresample/utils.py | 80 +++++---
pyresample/version.py | 2 +-
10 files changed, 709 insertions(+), 193 deletions(-)
diff --git a/.bumpversion.cfg b/.bumpversion.cfg
index eed2dcc..9232afd 100644
--- a/.bumpversion.cfg
+++ b/.bumpversion.cfg
@@ -1,5 +1,5 @@
[bumpversion]
-current_version = 1.6.0
+current_version = 1.6.1
commit = True
tag = True
diff --git a/changelog.rst b/changelog.rst
index cd16cc8..820a32d 100644
--- a/changelog.rst
+++ b/changelog.rst
@@ -2,6 +2,56 @@ Changelog
=========
+v1.6.1 (2017-09-18)
+-------------------
+- update changelog. [Martin Raspaud]
+- Bump version: 1.6.0 → 1.6.1. [Martin Raspaud]
+- Merge pull request #60 from pytroll/feature-dynamic-area. [David
+ Hoese]
+
+ Add support for dynamic areas
+- Merge branch 'develop' into feature-dynamic-area. [Martin Raspaud]
+- Apply assert_allclose to proj dicts for tests. [Martin Raspaud]
+- Fix some style issues. [Martin Raspaud]
+- Set DynamicArea proj to `omerc` by default. [Martin Raspaud]
+- Implement proposed changes in PR review. [Martin Raspaud]
+- Use numpy's assert almost equal for area_extent comparisons. [Martin
+ Raspaud]
+- Document the DynamicArea class. [Martin Raspaud]
+- Fix optimal projection computation tests. [Martin Raspaud]
+- Pep8 cleanup. [Martin Raspaud]
+- Valid index computation optimization. [Martin Raspaud]
+- Change bb computation api to use the whole proj_dict. [Martin Raspaud]
+- Fix unittests for updated omerc computations. [Martin Raspaud]
+- Use other azimuth direction for omerc. [Martin Raspaud]
+- Flip x and y size in omerc projection. [Martin Raspaud]
+- Bugfix typo. [Martin Raspaud]
+- Allow lons and lats to be any array in bb computation. [Martin
+ Raspaud]
+- Add SwathDefinition tests to the test suite. [Martin Raspaud]
+- Support bounding box area computation from SwathDefintion. [Martin
+ Raspaud]
+
+ This add support for computing a bounding box area from a swath definition that would fit optimally. The default projection is oblique mercator, with is optimal for locally received imager passes.
+- Add support for dynamic areas. [Martin Raspaud]
+- Merge pull request #70 from pytroll/feature-radius-parameters. [David
+ Hoese]
+
+ Add 'proj4_radius_parameters' to calculate 'a' and 'b' from ellps
+- Add tests for proj4_radius_parameters. [davidh-ssec]
+- Fix typo in function call in radius parameters. [davidh-ssec]
+- Add 'proj4_radius_parameters' to calculate 'a' and 'b' from ellps.
+ [davidh-ssec]
+- Merge pull request #68 from pytroll/feature-56. [Martin Raspaud]
+
+ Fix GridDefinition as permitted definition in preprocessing utils
+- Add more preprocessing tests. [davidh-ssec]
+- Fix preprocessing functions to use duck type on provided areas.
+ [davidh-ssec]
+- Fix GridDefinition as permitted definition in preprocessing utils.
+ [davidh-ssec]
+
+
v1.6.0 (2017-09-12)
-------------------
- update changelog. [Martin Raspaud]
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index c0f02ca..31d84ee 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -24,12 +24,16 @@
import warnings
from collections import OrderedDict
+from logging import getLogger
import numpy as np
import yaml
+from pyproj import Geod, Proj
from pyresample import _spatial_mp, utils
+logger = getLogger(__name__)
+
class DimensionError(Exception):
pass
@@ -59,6 +63,8 @@ class BaseDefinition(object):
if type(lons) != type(lats):
raise TypeError('lons and lats must be of same type')
elif lons is not None:
+ lons = np.asanyarray(lons)
+ lats = np.asanyarray(lats)
if lons.shape != lats.shape:
raise ValueError('lons and lats must have same shape')
@@ -351,6 +357,8 @@ class CoordinateDefinition(BaseDefinition):
"""Base class for geometry definitions defined by lons and lats only"""
def __init__(self, lons, lats, nprocs=1):
+ lons = np.asanyarray(lons)
+ lats = np.asanyarray(lats)
super(CoordinateDefinition, self).__init__(lons, lats, nprocs)
if lons.shape == lats.shape and lons.dtype == lats.dtype:
self.shape = lons.shape
@@ -449,12 +457,175 @@ class SwathDefinition(CoordinateDefinition):
"""
def __init__(self, lons, lats, nprocs=1):
+ lons = np.asanyarray(lons)
+ lats = np.asanyarray(lats)
super(SwathDefinition, self).__init__(lons, lats, nprocs)
if lons.shape != lats.shape:
raise ValueError('lon and lat arrays must have same shape')
elif lons.ndim > 2:
raise ValueError('Only 1 and 2 dimensional swaths are allowed')
+ def _compute_omerc_parameters(self, ellipsoid):
+ """Compute the oblique mercator projection bouding box parameters."""
+ lines, cols = self.lons.shape
+ lon1, lon2 = np.asanyarray(self.lons[[0, -1], int(cols / 2)])
+ lat1, lat, lat2 = np.asanyarray(
+ self.lats[[0, int(lines / 2), -1], int(cols / 2)])
+
+ proj_dict2points = {'proj': 'omerc', 'lat_0': lat, 'ellps': ellipsoid,
+ 'lat_1': lat1, 'lon_1': lon1,
+ 'lat_2': lat2, 'lon_2': lon2}
+
+ lonc, lat0 = Proj(**proj_dict2points)(0, 0, inverse=True)
+ az1, az2, dist = Geod(**proj_dict2points).inv(lonc, lat0, lon1, lat1)
+ del az2, dist
+ return {'proj': 'omerc', 'alpha': float(az1),
+ 'lat_0': float(lat0), 'lonc': float(lonc),
+ 'no_rot': True, 'ellps': ellipsoid}
+
+ def get_edge_lonlats(self):
+ """Get the concatenated boundary of the current swath."""
+ lons, lats = self.get_boundary_lonlats()
+ blons = np.ma.concatenate([lons.side1, lons.side2,
+ lons.side3, lons.side4])
+ blats = np.ma.concatenate([lats.side1, lats.side2,
+ lats.side3, lats.side4])
+ return blons, blats
+
+ def compute_bb_proj_params(self, proj_dict):
+ projection = proj_dict['proj']
+ ellipsoid = proj_dict.get('ellps', 'WGS84')
+ if projection == 'omerc':
+ return self._compute_omerc_parameters(ellipsoid)
+ else:
+ raise NotImplementedError('Only omerc supported for now.')
+
+ def compute_optimal_bb_area(self, proj_dict=None):
+ """Compute the "best" bounding box area for this swath with `proj_dict`.
+
+ By default, the projection is Oblique Mercator (`omerc` in proj.4), in
+ 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.
+ """
+ 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
+ x_size = int(cols * 1.1)
+ y_size = int(lines * 1.1)
+
+ proj_dict = self.compute_bb_proj_params(proj_dict)
+
+ if projection == 'omerc':
+ x_size, y_size = y_size, x_size
+
+ area = DynamicAreaDefinition(area_id, description, proj_dict)
+ lons, lats = self.get_edge_lonlats()
+ return area.freeze((lons, lats), size=(x_size, y_size))
+
+
+class DynamicAreaDefinition(object):
+ """An AreaDefintion containing just a subset of the needed parameters.
+
+ The purpose of this class is to be able to adapt the area extent and size of
+ the area to a given set of longitudes and latitudes, such that e.g. polar
+ satellite granules can be resampled optimaly to a give projection."""
+
+ def __init__(self, area_id=None, description=None, proj_dict=None,
+ x_size=None, y_size=None, area_extent=None,
+ optimize_projection=False):
+ """Initialize the DynamicAreaDefinition.
+
+ area_id:
+ The name of the area.
+ description:
+ The description of the area.
+ proj_dict:
+ The dictionary of projection parameters. Doesn't have to be complete.
+ x_size, y_size:
+ The size of the resulting area.
+ area_extent:
+ The area extent of the area.
+ optimize_projection:
+ Whether the projection parameters have to be optimized.
+ """
+ self.area_id = area_id
+ self.description = description
+ self.proj_dict = proj_dict
+ self.x_size = x_size
+ self.y_size = y_size
+ self.area_extent = area_extent
+ self.optimize_projection = optimize_projection
+
+ def compute_domain(self, corners, resolution=None, size=None):
+ """Compute size and area_extent from corners and [size or resolution]
+ info."""
+ if resolution is not None and size is not None:
+ raise ValueError("Both resolution and size can't be provided.")
+
+ if size:
+ x_size, y_size = size
+ x_resolution = (corners[2] - corners[0]) * 1.0 / (x_size - 1)
+ y_resolution = (corners[3] - corners[1]) * 1.0 / (y_size - 1)
+
+ if resolution:
+ try:
+ x_resolution, y_resolution = resolution
+ except TypeError:
+ x_resolution = y_resolution = resolution
+ x_size = int(np.rint((corners[2] - corners[0]) * 1.0 /
+ x_resolution + 1))
+ y_size = int(np.rint((corners[3] - corners[1]) * 1.0 /
+ y_resolution + 1))
+
+ area_extent = (corners[0] - x_resolution / 2,
+ corners[1] - y_resolution / 2,
+ corners[2] + x_resolution / 2,
+ corners[3] + y_resolution / 2)
+ return area_extent, x_size, y_size
+
+ def freeze(self, lonslats=None,
+ resolution=None, size=None,
+ proj_info=None):
+ """Create an AreaDefintion from this area with help of some extra info.
+
+ lonlats:
+ the geographical coordinates to contain in the resulting area.
+ resolution:
+ the resolution of the resulting area.
+ size:
+ the size of the resulting area.
+ proj_info:
+ complementing parameters to the projection info.
+
+ Resolution and Size parameters are ignored if the instance is created
+ with the `optimize_projection` flag set to True.
+ """
+ if proj_info is not None:
+ self.proj_dict.update(proj_info)
+
+ if self.optimize_projection:
+ return lonslats.compute_optimal_bb_area(self.proj_dict)
+
+ if not self.area_extent or not self.x_size or not self.y_size:
+ proj4 = Proj(**self.proj_dict)
+ try:
+ lons, lats = lonslats
+ except (TypeError, ValueError):
+ lons, lats = lonslats.get_lonlats()
+ xarr, yarr = proj4(np.asarray(lons), np.asarray(lats))
+ corners = [np.min(xarr), np.min(yarr), np.max(xarr), np.max(yarr)]
+
+ domain = self.compute_domain(corners, resolution, size)
+ self.area_extent, self.x_size, self.y_size = domain
+
+ return AreaDefinition(self.area_id, self.description, None,
+ self.proj_dict, self.x_size, self.y_size,
+ self.area_extent)
+
class AreaDefinition(BaseDefinition):
diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py
index 4be4322..e8cdcd2 100644
--- a/pyresample/kd_tree.py
+++ b/pyresample/kd_tree.py
@@ -21,13 +21,13 @@ supported"""
from __future__ import absolute_import
+import sys
import types
import warnings
-import sys
import numpy as np
-from pyresample import geometry, data_reduce, _spatial_mp
+from pyresample import _spatial_mp, data_reduce, geometry
if sys.version < '3':
range = xrange
@@ -66,20 +66,20 @@ def resample_nearest(source_geo_def, data, target_geo_def,
----------
source_geo_def : object
Geometry definition of source
- data : numpy array
+ data : numpy array
1d array of single channel data points or
(source_size, k) array of k channels of datapoints
target_geo_def : object
Geometry definition of target
- radius_of_influence : float
+ radius_of_influence : float
Cut off distance in meters
epsilon : float, optional
Allowed uncertainty in meters. Increasing uncertainty
reduces execution time
fill_value : int or None, optional
Set undetermined pixels to this value.
- If fill_value is None a masked array is returned
- with undetermined pixels masked
+ If fill_value is None a masked array is returned
+ with undetermined pixels masked
reduce_data : bool, optional
Perform initial coarse reduction of source dataset in order
to reduce execution time
@@ -91,7 +91,7 @@ def resample_nearest(source_geo_def, data, target_geo_def,
Returns
-------
- data : numpy array
+ data : numpy array
Source data resampled to target geometry
"""
@@ -110,26 +110,26 @@ def resample_gauss(source_geo_def, data, target_geo_def,
----------
source_geo_def : object
Geometry definition of source
- data : numpy array
+ data : numpy array
Array of single channel data points or
(source_geo_def.shape, k) array of k channels of datapoints
target_geo_def : object
Geometry definition of target
- radius_of_influence : float
+ radius_of_influence : float
Cut off distance in meters
- sigmas : list of floats or float
- List of sigmas to use for the gauss weighting of each
+ sigmas : list of floats or float
+ List of sigmas to use for the gauss weighting of each
channel 1 to k, w_k = exp(-dist^2/sigma_k^2).
If only one channel is resampled sigmas is a single float value.
- neighbours : int, optional
+ neighbours : int, optional
The number of neigbours to consider for each grid point
epsilon : float, optional
Allowed uncertainty in meters. Increasing uncertainty
reduces execution time
- fill_value : {int, None}, optional
+ fill_value : {int, None}, optional
Set undetermined pixels to this value.
- If fill_value is None a masked array is returned
- with undetermined pixels masked
+ If fill_value is None a masked array is returned
+ with undetermined pixels masked
reduce_data : bool, optional
Perform initial coarse reduction of source dataset in order
to reduce execution time
@@ -148,7 +148,7 @@ def resample_gauss(source_geo_def, data, target_geo_def,
data, stddev, counts : numpy array, numpy array, numpy array (if with_uncert == True)
Source data resampled to target geometry.
Weighted standard devaition for all pixels having more than one source value
- Counts of number of source values used in weighting per pixel
+ Counts of number of source values used in weighting per pixel
"""
def gauss(sigma):
@@ -190,27 +190,27 @@ def resample_custom(source_geo_def, data, target_geo_def,
----------
source_geo_def : object
Geometry definition of source
- data : numpy array
+ data : numpy array
Array of single channel data points or
(source_geo_def.shape, k) array of k channels of datapoints
target_geo_def : object
Geometry definition of target
- radius_of_influence : float
+ radius_of_influence : float
Cut off distance in meters
- weight_funcs : list of function objects or function object
- List of weight functions f(dist) to use for the weighting
+ weight_funcs : list of function objects or function object
+ List of weight functions f(dist) to use for the weighting
of each channel 1 to k.
If only one channel is resampled weight_funcs is
a single function object.
- neighbours : int, optional
+ neighbours : int, optional
The number of neigbours to consider for each grid point
epsilon : float, optional
Allowed uncertainty in meters. Increasing uncertainty
reduces execution time
- fill_value : {int, None}, optional
+ fill_value : {int, None}, optional
Set undetermined pixels to this value.
- If fill_value is None a masked array is returned
- with undetermined pixels masked
+ If fill_value is None a masked array is returned
+ with undetermined pixels masked
reduce_data : bool, optional
Perform initial coarse reduction of source dataset in order
to reduce execution time
@@ -282,9 +282,9 @@ def get_neighbour_info(source_geo_def, target_geo_def, radius_of_influence,
Geometry definition of source
target_geo_def : object
Geometry definition of target
- radius_of_influence : float
+ radius_of_influence : float
Cut off distance in meters
- neighbours : int, optional
+ neighbours : int, optional
The number of neigbours to consider for each grid point
epsilon : float, optional
Allowed uncertainty in meters. Increasing uncertainty
@@ -300,7 +300,7 @@ def get_neighbour_info(source_geo_def, target_geo_def, radius_of_influence,
Returns
-------
- (valid_input_index, valid_output_index,
+ (valid_input_index, valid_output_index,
index_array, distance_array) : tuple of numpy arrays
Neighbour resampling info
"""
@@ -390,8 +390,8 @@ def _get_valid_input_index(source_geo_def, target_geo_def, reduce_data,
"""Find indices of reduced inputput data"""
source_lons, source_lats = source_geo_def.get_lonlats(nprocs=nprocs)
- source_lons = source_lons.ravel()
- source_lats = source_lats.ravel()
+ source_lons = np.asanyarray(source_lons).ravel()
+ source_lats = np.asanyarray(source_lats).ravel()
if source_lons.size == 0 or source_lats.size == 0:
raise ValueError('Cannot resample empty data set')
@@ -400,9 +400,8 @@ def _get_valid_input_index(source_geo_def, target_geo_def, reduce_data,
raise ValueError('Mismatch between lons and lats')
# Remove illegal values
- valid_data = ((source_lons >= -180) & (source_lons <= 180) &
- (source_lats <= 90) & (source_lats >= -90))
- valid_input_index = np.ones(source_geo_def.size, dtype=np.bool)
+ valid_input_index = ((source_lons >= -180) & (source_lons <= 180) &
+ (source_lats <= 90) & (source_lats >= -90))
if reduce_data:
# Reduce dataset
@@ -415,16 +414,15 @@ def _get_valid_input_index(source_geo_def, target_geo_def, reduce_data,
geometry.AreaDefinition))):
# Resampling from swath to grid or from grid to grid
lonlat_boundary = target_geo_def.get_boundary_lonlats()
- valid_input_index = \
+
+ # Combine reduced and legal values
+ valid_input_index &= \
data_reduce.get_valid_index_from_lonlat_boundaries(
lonlat_boundary[0],
lonlat_boundary[1],
source_lons, source_lats,
radius_of_influence)
- # Combine reduced and legal values
- valid_input_index = (valid_data & valid_input_index)
-
if(isinstance(valid_input_index, np.ma.core.MaskedArray)):
# Make sure valid_input_index is not a masked array
valid_input_index = valid_input_index.filled(False)
@@ -469,7 +467,7 @@ def _create_resample_kdtree(source_lons, source_lats, valid_input_index, nprocs=
"""
if not isinstance(source_geo_def, geometry.BaseDefinition):
raise TypeError('source_geo_def must be of geometry type')
-
+
#Get reduced cartesian coordinates and flatten them
source_cartesian_coords = source_geo_def.get_cartesian_coords(nprocs=nprocs)
input_coords = geometry._flatten_cartesian_coords(source_cartesian_coords)
@@ -599,20 +597,20 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
distance_array : numpy array, optional
distance_array from get_neighbour_info
Not needed for 'nn' resample type
- weight_funcs : list of function objects or function object, optional
- List of weight functions f(dist) to use for the weighting
+ weight_funcs : list of function objects or function object, optional
+ List of weight functions f(dist) to use for the weighting
of each channel 1 to k.
If only one channel is resampled weight_funcs is
a single function object.
Must be supplied when using 'custom' resample type
fill_value : int or None, optional
Set undetermined pixels to this value.
- If fill_value is None a masked array is returned
+ If fill_value is None a masked array is returned
with undetermined pixels masked
Returns
-------
- result : numpy array
+ result : numpy array
Source data resampled to target geometry
"""
diff --git a/pyresample/test/test_ewa_ll2cr.py b/pyresample/test/test_ewa_ll2cr.py
index 99caaaa..018b0e3 100644
--- a/pyresample/test/test_ewa_ll2cr.py
+++ b/pyresample/test/test_ewa_ll2cr.py
@@ -24,6 +24,7 @@
import sys
import logging
import numpy as np
+from pyresample.test.utils import create_test_longitude, create_test_latitude
if sys.version_info < (2, 7):
import unittest2 as unittest
else:
@@ -32,28 +33,6 @@ else:
LOG = logging.getLogger(__name__)
-def create_test_longitude(start, stop, shape, twist_factor=0.0, dtype=np.float32):
- if start > stop:
- stop += 360.0
-
- lon_row = np.linspace(start, stop, num=shape[1]).astype(dtype)
- twist_array = np.arange(shape[0]).reshape((shape[0], 1)) * twist_factor
- lon_array = np.repeat([lon_row], shape[0], axis=0)
- lon_array += twist_array
-
- if stop > 360.0:
- lon_array[lon_array > 360.0] -= 360
- return lon_array
-
-
-def create_test_latitude(start, stop, shape, twist_factor=0.0, dtype=np.float32):
- lat_col = np.linspace(start, stop, num=shape[0]).astype(dtype).reshape((shape[0], 1))
- twist_array = np.arange(shape[1]) * twist_factor
- lat_array = np.repeat(lat_col, shape[1], axis=1)
- lat_array += twist_array
- return lat_array
-
-
dynamic_wgs84 = {
"grid_name": "test_wgs84_fit",
"origin_x": None,
diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py
index 164d32c..73a3050 100644
--- a/pyresample/test/test_geometry.py
+++ b/pyresample/test/test_geometry.py
@@ -165,39 +165,6 @@ class Test(unittest.TestCase):
"BaseDefinition did not maintain dtype of longitudes (in:%s out:%s)" %
(lons2_ints.dtype, lons.dtype,))
- def test_swath(self):
- lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
- lats1 = np.fromfunction(
- lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
-
- swath_def = geometry.SwathDefinition(lons1, lats1)
-
- lons2, lats2 = swath_def.get_lonlats()
-
- self.assertFalse(id(lons1) != id(lons2) or id(lats1) != id(lats2),
- msg='Caching of swath coordinates failed')
-
- def test_swath_wrap(self):
- lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
- lats1 = np.fromfunction(
- lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
-
- lons1 += 180.
- with catch_warnings() as w1:
- swath_def = geometry.BaseDefinition(lons1, lats1)
- self.assertFalse(
- len(w1) != 1, 'Failed to trigger a warning on longitude wrapping')
- self.assertFalse(('-180:+180' not in str(w1[0].message)),
- 'Failed to trigger correct warning about longitude wrapping')
-
- lons2, lats2 = swath_def.get_lonlats()
-
- self.assertTrue(id(lons1) != id(lons2),
- msg='Caching of swath coordinates failed with longitude wrapping')
-
- self.assertTrue(lons2.min() > -180 and lons2.max() < 180,
- 'Wrapping of longitudes failed for SwathDefinition')
-
def test_area_equal(self):
area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
{'a': '6378144.0',
@@ -260,23 +227,6 @@ class Test(unittest.TestCase):
self.assertFalse(
area_def == msg_area, 'area_defs are not expected to be equal')
- def test_swath_equal(self):
- lons = np.array([1.2, 1.3, 1.4, 1.5])
- lats = np.array([65.9, 65.86, 65.82, 65.78])
- swath_def = geometry.SwathDefinition(lons, lats)
- swath_def2 = geometry.SwathDefinition(lons, lats)
- self.assertFalse(
- swath_def != swath_def2, 'swath_defs are not equal as expected')
-
- def test_swath_not_equal(self):
- lats1 = np.array([65.9, 65.86, 65.82, 65.78])
- lons = np.array([1.2, 1.3, 1.4, 1.5])
- lats2 = np.array([65.91, 65.85, 65.80, 65.75])
- swath_def = geometry.SwathDefinition(lons, lats1)
- swath_def2 = geometry.SwathDefinition(lons, lats2)
- self.assertFalse(
- swath_def == swath_def2, 'swath_defs are not expected to be equal')
-
def test_swath_equal_area(self):
area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
{'a': '6378144.0',
@@ -353,60 +303,6 @@ class Test(unittest.TestCase):
self.assertFalse(
area_def == swath_def, "swath_def and area_def should be different")
- def test_concat_1d(self):
- lons1 = np.array([1, 2, 3])
- lats1 = np.array([1, 2, 3])
- lons2 = np.array([4, 5, 6])
- lats2 = np.array([4, 5, 6])
- swath_def1 = geometry.SwathDefinition(lons1, lats1)
- swath_def2 = geometry.SwathDefinition(lons2, lats2)
- swath_def_concat = swath_def1.concatenate(swath_def2)
- expected = np.array([1, 2, 3, 4, 5, 6])
- self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and
- np.array_equal(swath_def_concat.lons, expected),
- 'Failed to concatenate 1D swaths')
-
- def test_concat_2d(self):
- lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
- lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
- lons2 = np.array([[4, 5, 6], [6, 7, 8]])
- lats2 = np.array([[4, 5, 6], [6, 7, 8]])
- swath_def1 = geometry.SwathDefinition(lons1, lats1)
- swath_def2 = geometry.SwathDefinition(lons2, lats2)
- swath_def_concat = swath_def1.concatenate(swath_def2)
- expected = np.array(
- [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]])
- self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and
- np.array_equal(swath_def_concat.lons, expected),
- 'Failed to concatenate 2D swaths')
-
- def test_append_1d(self):
- lons1 = np.array([1, 2, 3])
- lats1 = np.array([1, 2, 3])
- lons2 = np.array([4, 5, 6])
- lats2 = np.array([4, 5, 6])
- swath_def1 = geometry.SwathDefinition(lons1, lats1)
- swath_def2 = geometry.SwathDefinition(lons2, lats2)
- swath_def1.append(swath_def2)
- expected = np.array([1, 2, 3, 4, 5, 6])
- self.assertTrue(np.array_equal(swath_def1.lons, expected) and
- np.array_equal(swath_def1.lons, expected),
- 'Failed to append 1D swaths')
-
- def test_append_2d(self):
- lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
- lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
- lons2 = np.array([[4, 5, 6], [6, 7, 8]])
- lats2 = np.array([[4, 5, 6], [6, 7, 8]])
- swath_def1 = geometry.SwathDefinition(lons1, lats1)
- swath_def2 = geometry.SwathDefinition(lons2, lats2)
- swath_def1.append(swath_def2)
- expected = np.array(
- [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]])
- self.assertTrue(np.array_equal(swath_def1.lons, expected) and
- np.array_equal(swath_def1.lons, expected),
- 'Failed to append 2D swaths')
-
def test_grid_filter_valid(self):
lons = np.array([-170, -30, 30, 170])
lats = np.array([20, -40, 50, -80])
@@ -706,6 +602,205 @@ class Test(unittest.TestCase):
self.assertTrue((y__.data == y_expects).all())
+def assert_np_dict_allclose(dict1, dict2):
+
+ assert set(dict1.keys()) == set(dict2.keys())
+ for key, val in dict1.items():
+ try:
+ np.testing.assert_allclose(val, dict2[key])
+ except TypeError:
+ assert(val == dict2[key])
+
+
+class TestSwathDefinition(unittest.TestCase):
+ """Test the SwathDefinition."""
+
+ def test_swath(self):
+ lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
+ lats1 = np.fromfunction(
+ lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
+
+ swath_def = geometry.SwathDefinition(lons1, lats1)
+
+ lons2, lats2 = swath_def.get_lonlats()
+
+ self.assertFalse(id(lons1) != id(lons2) or id(lats1) != id(lats2),
+ msg='Caching of swath coordinates failed')
+
+ def test_swath_wrap(self):
+ lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100))
+ lats1 = np.fromfunction(
+ lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100))
+
+ lons1 += 180.
+ with catch_warnings() as w1:
+ swath_def = geometry.BaseDefinition(lons1, lats1)
+ self.assertFalse(
+ len(w1) != 1, 'Failed to trigger a warning on longitude wrapping')
+ self.assertFalse(('-180:+180' not in str(w1[0].message)),
+ 'Failed to trigger correct warning about longitude wrapping')
+
+ lons2, lats2 = swath_def.get_lonlats()
+
+ self.assertTrue(id(lons1) != id(lons2),
+ msg='Caching of swath coordinates failed with longitude wrapping')
+
+ self.assertTrue(lons2.min() > -180 and lons2.max() < 180,
+ 'Wrapping of longitudes failed for SwathDefinition')
+
+ def test_concat_1d(self):
+ lons1 = np.array([1, 2, 3])
+ lats1 = np.array([1, 2, 3])
+ lons2 = np.array([4, 5, 6])
+ lats2 = np.array([4, 5, 6])
+ swath_def1 = geometry.SwathDefinition(lons1, lats1)
+ swath_def2 = geometry.SwathDefinition(lons2, lats2)
+ swath_def_concat = swath_def1.concatenate(swath_def2)
+ expected = np.array([1, 2, 3, 4, 5, 6])
+ self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and
+ np.array_equal(swath_def_concat.lons, expected),
+ 'Failed to concatenate 1D swaths')
+
+ def test_concat_2d(self):
+ lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
+ lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
+ lons2 = np.array([[4, 5, 6], [6, 7, 8]])
+ lats2 = np.array([[4, 5, 6], [6, 7, 8]])
+ swath_def1 = geometry.SwathDefinition(lons1, lats1)
+ swath_def2 = geometry.SwathDefinition(lons2, lats2)
+ swath_def_concat = swath_def1.concatenate(swath_def2)
+ expected = np.array(
+ [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]])
+ self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and
+ np.array_equal(swath_def_concat.lons, expected),
+ 'Failed to concatenate 2D swaths')
+
+ def test_append_1d(self):
+ lons1 = np.array([1, 2, 3])
+ lats1 = np.array([1, 2, 3])
+ lons2 = np.array([4, 5, 6])
+ lats2 = np.array([4, 5, 6])
+ swath_def1 = geometry.SwathDefinition(lons1, lats1)
+ swath_def2 = geometry.SwathDefinition(lons2, lats2)
+ swath_def1.append(swath_def2)
+ expected = np.array([1, 2, 3, 4, 5, 6])
+ self.assertTrue(np.array_equal(swath_def1.lons, expected) and
+ np.array_equal(swath_def1.lons, expected),
+ 'Failed to append 1D swaths')
+
+ def test_append_2d(self):
+ lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
+ lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]])
+ lons2 = np.array([[4, 5, 6], [6, 7, 8]])
+ lats2 = np.array([[4, 5, 6], [6, 7, 8]])
+ swath_def1 = geometry.SwathDefinition(lons1, lats1)
+ swath_def2 = geometry.SwathDefinition(lons2, lats2)
+ swath_def1.append(swath_def2)
+ expected = np.array(
+ [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]])
+ self.assertTrue(np.array_equal(swath_def1.lons, expected) and
+ np.array_equal(swath_def1.lons, expected),
+ 'Failed to append 2D swaths')
+
+ def test_swath_equal(self):
+ """Test swath equality."""
+ lons = np.array([1.2, 1.3, 1.4, 1.5])
+ lats = np.array([65.9, 65.86, 65.82, 65.78])
+ swath_def = geometry.SwathDefinition(lons, lats)
+ swath_def2 = geometry.SwathDefinition(lons, lats)
+ self.assertFalse(
+ swath_def != swath_def2, 'swath_defs are not equal as expected')
+
+ def test_swath_not_equal(self):
+ """Test swath inequality."""
+ lats1 = np.array([65.9, 65.86, 65.82, 65.78])
+ lons = np.array([1.2, 1.3, 1.4, 1.5])
+ lats2 = np.array([65.91, 65.85, 65.80, 65.75])
+ swath_def = geometry.SwathDefinition(lons, lats1)
+ swath_def2 = geometry.SwathDefinition(lons, lats2)
+ self.assertFalse(
+ swath_def == swath_def2, 'swath_defs are not expected to be equal')
+
+ def test_compute_omerc_params(self):
+ """Test omerc parameters computation."""
+ lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+ [80.84000396728516, 60.74200439453125, 34.08500289916992],
+ [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+
+ lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+ [79.11000061035156, 7.284000396728516, -5.107000350952148],
+ [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+
+ area = geometry.SwathDefinition(lons, lats)
+ proj_dict = {'no_rot': True, 'lonc': 5.340645620216994,
+ 'ellps': 'WGS84', 'proj': 'omerc',
+ 'alpha': 19.022450179020247, 'lat_0': 60.7420043944989}
+ assert_np_dict_allclose(area._compute_omerc_parameters('WGS84'),
+ proj_dict)
+
+ def test_get_edge_lonlats(self):
+ """Test the `get_edge_lonlats` functionality."""
+ lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+ [80.84000396728516, 60.74200439453125, 34.08500289916992],
+ [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+
+ lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+ [79.11000061035156, 7.284000396728516, -5.107000350952148],
+ [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+
+ area = geometry.SwathDefinition(lons, lats)
+ lons, lats = area.get_edge_lonlats()
+
+ np.testing.assert_allclose(lons, [-90.67900085, 79.11000061, 81.26400757,
+ 81.26400757, 29.67200089, 10.26000023,
+ 10.26000023, -5.10700035, -21.52500153,
+ -21.52500153, -21.56500053, -90.67900085])
+ np.testing.assert_allclose(lats, [85.23900604, 80.84000397, 67.07600403,
+ 67.07600403, 54.14700317, 30.54700089,
+ 30.54700089, 34.0850029, 35.58000183,
+ 35.58000183, 62.25600433, 85.23900604])
+
+ lats = np.array([[80., 80., 80.],
+ [80., 90., 80],
+ [80., 80., 80.]]).T
+
+ lons = np.array([[-45., 0., 45.],
+ [-90, 0., 90.],
+ [-135., 180., 135.]]).T
+
+ area = geometry.SwathDefinition(lons, lats)
+ lons, lats = area.get_edge_lonlats()
+
+ np.testing.assert_allclose(lons, [-45., -90., -135., -135., -180., 135.,
+ 135., 90., 45., 45., 0., -45.])
+ np.testing.assert_allclose(lats, [80., 80., 80., 80., 80., 80., 80.,
+ 80., 80., 80., 80., 80.])
+
+ def test_compute_optimal_bb(self):
+ """Test computing the bb area."""
+ lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469],
+ [80.84000396728516, 60.74200439453125, 34.08500289916992],
+ [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T
+
+ lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906],
+ [79.11000061035156, 7.284000396728516, -5.107000350952148],
+ [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T
+
+ area = geometry.SwathDefinition(lons, lats)
+
+ res = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'})
+
+ np.testing.assert_allclose(res.area_extent, (2286629.731529,
+ -2359693.817959,
+ 11729881.856072,
+ 2437001.523925))
+ proj_dict = {'no_rot': True, 'lonc': 5.340645620216994,
+ 'ellps': 'WGS84', 'proj': 'omerc',
+ 'alpha': 19.022450179020247, 'lat_0': 60.7420043944989}
+ assert_np_dict_allclose(res.proj_dict, proj_dict)
+ self.assertEqual(res.shape, (3, 3))
+
+
class TestStackedAreaDefinition(unittest.TestCase):
"""Test the StackedAreaDefition."""
@@ -864,6 +959,86 @@ class TestStackedAreaDefinition(unittest.TestCase):
area_extent)
+class TestDynamicAreaDefinition(unittest.TestCase):
+ """Test the DynamicAreaDefinition class."""
+
+ def test_freeze(self):
+ """Test freezing the area."""
+ area = geometry.DynamicAreaDefinition('test_area', 'A test area',
+ {'proj': 'laea'})
+ lons = [10, 10, 22, 22]
+ lats = [50, 66, 66, 50]
+ result = area.freeze((lons, lats),
+ resolution=3000,
+ proj_info={'lon0': 16, 'lat0': 58})
+
+ np.testing.assert_allclose(result.area_extent, (538546.7274949469,
+ 5380808.879250369,
+ 1724415.6519203288,
+ 6998895.701001488))
+ self.assertEqual(result.proj_dict['lon0'], 16)
+ self.assertEqual(result.proj_dict['lat0'], 58)
+ self.assertEqual(result.x_size, 395)
+ self.assertEqual(result.y_size, 539)
+
+ def test_freeze_with_bb(self):
+ """Test freezing the area with bounding box computation."""
+ # area = geometry.DynamicAreaDefinition('test_area', 'A test area',
+ # {'proj': 'omerc'},
+ # optimize_projection=False)
+ # lons = [[10, 12.1, 14.2, 16.3],
+ # [10, 12, 14, 16],
+ # [10, 11.9, 13.8, 15.7]]
+ # lats = [[66, 67, 68, 69.],
+ # [58, 59, 60, 61],
+ # [50, 51, 52, 53]]
+ # sdef = geometry.SwathDefinition(lons, lats)
+ # result = area.freeze(sdef,
+ # resolution=1000)
+ # self.assertTupleEqual(result.area_extent, (5578795.1654752363,
+ # -270848.61872542271,
+ # 7694893.3964453982,
+ # 126974.877141819))
+ # self.assertEqual(result.x_size, 2116)
+ # self.assertEqual(result.y_size, 398)
+
+ area = geometry.DynamicAreaDefinition('test_area', 'A test area',
+ {'proj': 'omerc'},
+ optimize_projection=True)
+ lons = [[10, 12.1, 14.2, 16.3],
+ [10, 12, 14, 16],
+ [10, 11.9, 13.8, 15.7]]
+ lats = [[66, 67, 68, 69.],
+ [58, 59, 60, 61],
+ [50, 51, 52, 53]]
+ sdef = geometry.SwathDefinition(lons, lats)
+ result = area.freeze(sdef,
+ resolution=1000)
+ np.testing.assert_allclose(result.area_extent, (5050520.6077326955,
+ -336485.86803662963,
+ 8223167.9541879389,
+ 192612.12645302597))
+ self.assertEqual(result.x_size, 3)
+ self.assertEqual(result.y_size, 4)
+
+ def test_compute_domain(self):
+ """Test computing size and area extent."""
+ area = geometry.DynamicAreaDefinition('test_area', 'A test area',
+ {'proj': 'laea'})
+ corners = [1, 1, 9, 9]
+ self.assertRaises(ValueError, area.compute_domain, corners, 1, 1)
+
+ area_extent, x_size, y_size = area.compute_domain(corners, size=(5, 5))
+ self.assertTupleEqual(area_extent, (0, 0, 10, 10))
+ self.assertEqual(x_size, 5)
+ self.assertEqual(y_size, 5)
+
+ area_extent, x_size, y_size = area.compute_domain(corners, resolution=2)
+ self.assertTupleEqual(area_extent, (0, 0, 10, 10))
+ self.assertEqual(x_size, 5)
+ self.assertEqual(y_size, 5)
+
+
def suite():
"""The test suite.
"""
@@ -871,6 +1046,8 @@ def suite():
mysuite = unittest.TestSuite()
mysuite.addTest(loader.loadTestsFromTestCase(Test))
mysuite.addTest(loader.loadTestsFromTestCase(TestStackedAreaDefinition))
+ mysuite.addTest(loader.loadTestsFromTestCase(TestDynamicAreaDefinition))
+ mysuite.addTest(loader.loadTestsFromTestCase(TestSwathDefinition))
return mysuite
diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py
index d79af05..7f0409e 100644
--- a/pyresample/test/test_utils.py
+++ b/pyresample/test/test_utils.py
@@ -3,7 +3,7 @@ import unittest
import numpy as np
-from pyresample import utils
+from pyresample.test.utils import create_test_longitude, create_test_latitude
def tmp(f):
@@ -14,6 +14,7 @@ def tmp(f):
class TestLegacyAreaParser(unittest.TestCase):
def test_area_parser_legacy(self):
"""Test legacy area parser."""
+ from pyresample import utils
ease_nh, ease_sh = utils.parse_area_file(os.path.join(os.path.dirname(__file__),
'test_files',
'areas.cfg'), 'ease_nh', 'ease_sh')
@@ -37,6 +38,7 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
self.assertEquals(ease_sh.__str__(), sh_str)
def test_load_area(self):
+ from pyresample import utils
ease_nh = utils.load_area(os.path.join(os.path.dirname(__file__),
'test_files',
'areas.cfg'), 'ease_nh')
@@ -50,6 +52,7 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
self.assertEquals(nh_str, ease_nh.__str__())
def test_not_found_exception(self):
+ from pyresample import utils
self.assertRaises(utils.AreaNotFound, utils.parse_area_file,
os.path.join(
os.path.dirname(__file__), 'test_files', 'areas.cfg'),
@@ -59,6 +62,7 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
class TestYAMLAreaParser(unittest.TestCase):
def test_area_parser_yaml(self):
"""Test YAML area parser."""
+ from pyresample import utils
ease_nh, ease_sh = utils.parse_area_file(os.path.join(os.path.dirname(__file__),
'test_files',
'areas.yaml'),
@@ -81,6 +85,7 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
self.assertEquals(ease_sh.__str__(), sh_str)
def test_multiple_file_content(self):
+ from pyresample import utils
area_list = ["""ease_sh:
description: Antarctic EASE grid
projection:
@@ -119,9 +124,62 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)"""
self.assertIn(results[1].area_id, ('ease_sh', 'ease_sh2'))
+class TestPreprocessing(unittest.TestCase):
+ def test_nearest_neighbor_area_area(self):
+ from pyresample import utils, geometry
+ proj_str = "+proj=lcc +datum=WGS84 +ellps=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs"
+ proj_dict = utils.proj4_str_to_dict(proj_str)
+ extents = [0, 0, 1000. * 5000, 1000. * 5000]
+ area_def = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS',
+ proj_dict, 400, 500, extents)
+
+ extents2 = [-1000, -1000, 1000. * 4000, 1000. * 4000]
+ area_def2 = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS',
+ proj_dict, 600, 700, extents2)
+ rows, cols = utils.generate_nearest_neighbour_linesample_arrays(area_def, area_def2, 12000.)
+
+ def test_nearest_neighbor_area_grid(self):
+ from pyresample import utils, geometry
+ lon_arr = create_test_longitude(-94.9, -90.0, (50, 100), dtype=np.float64)
+ lat_arr = create_test_latitude(25.1, 30.0, (50, 100), dtype=np.float64)
+ grid = geometry.GridDefinition(lons=lon_arr, lats=lat_arr)
+
+ proj_str = "+proj=lcc +datum=WGS84 +ellps=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs"
+ proj_dict = utils.proj4_str_to_dict(proj_str)
+ extents = [0, 0, 1000. * 5000, 1000. * 5000]
+ area_def = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS',
+ proj_dict, 400, 500, extents)
+ rows, cols = utils.generate_nearest_neighbour_linesample_arrays(area_def, grid, 12000.)
+
+ def test_nearest_neighbor_grid_area(self):
+ from pyresample import utils, geometry
+ proj_str = "+proj=lcc +datum=WGS84 +ellps=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs"
+ proj_dict = utils.proj4_str_to_dict(proj_str)
+ extents = [0, 0, 1000. * 2500., 1000. * 2000.]
+ area_def = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS',
+ proj_dict, 40, 50, extents)
+
+ lon_arr = create_test_longitude(-100.0, -60.0, (550, 500), dtype=np.float64)
+ lat_arr = create_test_latitude(20.0, 45.0, (550, 500), dtype=np.float64)
+ grid = geometry.GridDefinition(lons=lon_arr, lats=lat_arr)
+ rows, cols = utils.generate_nearest_neighbour_linesample_arrays(grid, area_def, 12000.)
+
+ def test_nearest_neighbor_grid_grid(self):
+ from pyresample import utils, geometry
+ lon_arr = create_test_longitude(-95.0, -85.0, (40, 50), dtype=np.float64)
+ lat_arr = create_test_latitude(25.0, 35.0, (40, 50), dtype=np.float64)
+ grid_dst = geometry.GridDefinition(lons=lon_arr, lats=lat_arr)
+
+ lon_arr = create_test_longitude(-100.0, -80.0, (400, 500), dtype=np.float64)
+ lat_arr = create_test_latitude(20.0, 40.0, (400, 500), dtype=np.float64)
+ grid = geometry.GridDefinition(lons=lon_arr, lats=lat_arr)
+ rows, cols = utils.generate_nearest_neighbour_linesample_arrays(grid, grid_dst, 12000.)
+
+
class TestMisc(unittest.TestCase):
def test_wrap_longitudes(self):
# test that we indeed wrap to [-180:+180[
+ from pyresample import utils
step = 60
lons = np.arange(-360, 360 + step, step)
self.assertTrue(
@@ -133,10 +191,36 @@ class TestMisc(unittest.TestCase):
def test_unicode_proj4_string(self):
"""Test that unicode is accepted for area creation.
"""
+ from pyresample import utils
utils.get_area_def(u"eurol", u"eurol", u"bla",
u'+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45',
1000, 1000, (-1000, -1000, 1000, 1000))
+ def test_proj4_radius_parameters_provided(self):
+ from pyresample import utils
+ a, b = utils.proj4_radius_parameters(
+ '+proj=stere +a=6378273 +b=6356889.44891',
+ )
+ np.testing.assert_almost_equal(a, 6378273)
+ np.testing.assert_almost_equal(b, 6356889.44891)
+
+ def test_proj4_radius_parameters_ellps(self):
+ from pyresample import utils
+ a, b = utils.proj4_radius_parameters(
+ '+proj=stere +ellps=WGS84',
+ )
+ np.testing.assert_almost_equal(a, 6378137.)
+ np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
+
+ def test_proj4_radius_parameters_default(self):
+ from pyresample import utils
+ a, b = utils.proj4_radius_parameters(
+ '+proj=lcc',
+ )
+ # WGS84
+ np.testing.assert_almost_equal(a, 6378137.)
+ np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
+
def suite():
"""The test suite.
@@ -145,6 +229,7 @@ def suite():
mysuite = unittest.TestSuite()
mysuite.addTest(loader.loadTestsFromTestCase(TestLegacyAreaParser))
mysuite.addTest(loader.loadTestsFromTestCase(TestYAMLAreaParser))
+ mysuite.addTest(loader.loadTestsFromTestCase(TestPreprocessing))
mysuite.addTest(loader.loadTestsFromTestCase(TestMisc))
return mysuite
diff --git a/pyresample/test/utils.py b/pyresample/test/utils.py
index 39033e4..1dbbcf9 100644
--- a/pyresample/test/utils.py
+++ b/pyresample/test/utils.py
@@ -25,6 +25,7 @@ import six
import types
import warnings
+import numpy as np
_deprecations_as_exceptions = False
_include_astropy_deprecations = False
@@ -151,3 +152,24 @@ class catch_warnings(warnings.catch_warnings):
treat_deprecations_as_exceptions()
+def create_test_longitude(start, stop, shape, twist_factor=0.0, dtype=np.float32):
+ if start > stop:
+ stop += 360.0
+
+ lon_row = np.linspace(start, stop, num=shape[1]).astype(dtype)
+ twist_array = np.arange(shape[0]).reshape((shape[0], 1)) * twist_factor
+ lon_array = np.repeat([lon_row], shape[0], axis=0)
+ lon_array += twist_array
+
+ if stop > 360.0:
+ lon_array[lon_array > 360.0] -= 360
+ return lon_array
+
+
+def create_test_latitude(start, stop, shape, twist_factor=0.0, dtype=np.float32):
+ lat_col = np.linspace(start, stop, num=shape[0]).astype(dtype).reshape((shape[0], 1))
+ twist_array = np.arange(shape[1]) * twist_factor
+ lat_array = np.repeat(lat_col, shape[1], axis=1)
+ lat_array += twist_array
+ return lat_array
+
diff --git a/pyresample/utils.py b/pyresample/utils.py
index 2bb70ed..3351b54 100644
--- a/pyresample/utils.py
+++ b/pyresample/utils.py
@@ -137,14 +137,27 @@ def _parse_yaml_area_file(area_file_name, *regions):
area_name, area_file_name))
description = params['description']
projection = params['projection']
- xsize = params['shape']['width']
- ysize = params['shape']['height']
- area_extent = (params['area_extent']['lower_left_xy'] +
- params['area_extent']['upper_right_xy'])
- res.append(pr.geometry.AreaDefinition(area_name, description,
- None, projection,
- xsize, ysize,
- area_extent))
+ optimize_projection = params.get('optimize_projection', False)
+ try:
+ xsize = params['shape']['width']
+ ysize = params['shape']['height']
+ except KeyError:
+ xsize, ysize = None, None
+ try:
+ area_extent = (params['area_extent']['lower_left_xy'] +
+ params['area_extent']['upper_right_xy'])
+ except KeyError:
+ area_extent = None
+ area = pr.geometry.DynamicAreaDefinition(area_name, description,
+ projection, xsize, ysize,
+ area_extent,
+ optimize_projection)
+ try:
+ area = area.freeze()
+ except (TypeError, AttributeError):
+ pass
+
+ res.append(area)
return res
@@ -278,9 +291,9 @@ def generate_quick_linesample_arrays(source_area_def, target_area_def, nprocs=1)
Parameters
-----------
source_area_def : object
- Source area definition as AreaDefinition object
+ Source area definition as geometry definition object
target_area_def : object
- Target area definition as AreaDefinition object
+ Target area definition as geometry definition object
nprocs : int, optional
Number of processor cores to be used
@@ -288,11 +301,6 @@ def generate_quick_linesample_arrays(source_area_def, target_area_def, nprocs=1)
-------
(row_indices, col_indices) : tuple of numpy arrays
"""
- if not (isinstance(source_area_def, pr.geometry.AreaDefinition) and
- isinstance(target_area_def, pr.geometry.AreaDefinition)):
- raise TypeError('source_area_def and target_area_def must be of type '
- 'geometry.AreaDefinition')
-
lons, lats = target_area_def.get_lonlats(nprocs)
source_pixel_y, source_pixel_x = pr.grid.get_linesample(lons, lats,
@@ -314,9 +322,9 @@ def generate_nearest_neighbour_linesample_arrays(source_area_def, target_area_de
Parameters
-----------
source_area_def : object
- Source area definition as AreaDefinition object
+ Source area definition as geometry definition object
target_area_def : object
- Target area definition as AreaDefinition object
+ Target area definition as geometry definition object
radius_of_influence : float
Cut off distance in meters
nprocs : int, optional
@@ -327,11 +335,6 @@ def generate_nearest_neighbour_linesample_arrays(source_area_def, target_area_de
(row_indices, col_indices) : tuple of numpy arrays
"""
- if not (isinstance(source_area_def, pr.geometry.AreaDefinition) and
- isinstance(target_area_def, pr.geometry.AreaDefinition)):
- raise TypeError('source_area_def and target_area_def must be of type '
- 'geometry.AreaDefinition')
-
valid_input_index, valid_output_index, index_array, distance_array = \
pr.kd_tree.get_neighbour_info(source_area_def,
target_area_def,
@@ -400,13 +403,44 @@ def _get_proj4_args(proj4_args):
def proj4_str_to_dict(proj4_str):
"""Convert PROJ.4 compatible string definition to dict
-
+
Note: Key only parameters will be assigned a value of `True`.
"""
pairs = (x.split('=', 1) for x in proj4_str.split(" "))
return dict((x[0], (x[1] if len(x) == 2 else True)) for x in pairs)
+def proj4_radius_parameters(proj4_dict):
+ """Calculate 'a' and 'b' radius parameters.
+
+ Arguments:
+ proj4_dict (str or dict): PROJ.4 parameters
+
+ Returns:
+ a (float), b (float): equatorial and polar radius
+ """
+ if isinstance(proj4_dict, str):
+ new_info = proj4_str_to_dict(proj4_dict)
+ else:
+ new_info = proj4_dict.copy()
+
+ # load information from PROJ.4 about the ellipsis if possible
+ if '+a' not in new_info or '+b' not in new_info:
+ import pyproj
+ ellps = pyproj.pj_ellps[new_info.get('+ellps', 'WGS84')]
+ new_info['+a'] = ellps['a']
+ if 'b' not in ellps and 'rf' in ellps:
+ new_info['+f'] = 1. / ellps['rf']
+ else:
+ new_info['+b'] = ellps['b']
+
+ if '+a' in new_info and '+f' in new_info and '+b' not in new_info:
+ # add a 'b' attribute back in if they used 'f' instead
+ new_info['+b'] = new_info['+a'] * (1 - new_info['+f'])
+
+ return float(new_info['+a']), float(new_info['+b'])
+
+
def _downcast_index_array(index_array, size):
"""Try to downcast array to uint16
"""
diff --git a/pyresample/version.py b/pyresample/version.py
index ee9fd28..6f605a3 100644
--- a/pyresample/version.py
+++ b/pyresample/version.py
@@ -15,4 +15,4 @@
# 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/>.
-__version__ = '1.6.0'
+__version__ = '1.6.1'
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-grass/pyresample.git
More information about the Pkg-grass-devel
mailing list