[pyresample] 01/05: New upstream version 1.7.0
Antonio Valentino
a_valentino-guest at moszumanska.debian.org
Sat Oct 14 07:57:20 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 89a964b091b7380226e5c23e545c5a7dcc762468
Author: Antonio Valentino <antonio.valentino at tiscali.it>
Date: Sat Oct 14 06:54:50 2017 +0000
New upstream version 1.7.0
---
.bumpversion.cfg | 2 +-
.travis.yml | 14 +-
appveyor.yml | 13 +-
changelog.rst | 97 ++++++++++++++
docs/source/swath.rst | 20 +--
pyresample/bilinear/__init__.py | 29 +++-
pyresample/geometry.py | 113 +++++++++++++---
pyresample/image.py | 184 +++++++++++++++++++++-----
pyresample/kd_tree.py | 277 +++++++++++++++++++++++++++++++++++++++
pyresample/test/test_bilinear.py | 6 +-
pyresample/test/test_image.py | 80 ++++++++---
pyresample/test/test_kd_tree.py | 12 +-
pyresample/test/test_utils.py | 8 ++
pyresample/utils.py | 88 +++++++++----
pyresample/version.py | 2 +-
15 files changed, 811 insertions(+), 134 deletions(-)
diff --git a/.bumpversion.cfg b/.bumpversion.cfg
index 9232afd..064551a 100644
--- a/.bumpversion.cfg
+++ b/.bumpversion.cfg
@@ -1,5 +1,5 @@
[bumpversion]
-current_version = 1.6.1
+current_version = 1.7.0
commit = True
tag = True
diff --git a/.travis.yml b/.travis.yml
index f433d0c..f7c6604 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,7 +1,6 @@
language: python
python:
- '2.7'
-- '3.3'
- '3.4'
- '3.5'
- '3.6'
@@ -11,21 +10,20 @@ before_install:
- sudo apt-get install libfreetype6-dev
- sudo apt-get install libgeos-dev
install:
-- if [[ $TRAVIS_PYTHON_VERSION == "2.7" ]]; then pip install "matplotlib>=1.5.0"; fi
+# matplotlib 2.1.0 has a bug that causes plotting masked arrays to fail
+# https://github.com/matplotlib/matplotlib/issues/9280
+- if [[ $TRAVIS_PYTHON_VERSION == "2.7" ]]; then pip install "matplotlib>=1.5.0,!=2.1.0"; fi
- if [[ $TRAVIS_PYTHON_VERSION == "2.7" ]]; then pip install "sphinx>=1.5.0"; fi
-- if [[ $TRAVIS_PYTHON_VERSION == "3.3" ]]; then pip install "matplotlib<1.5.0"; fi
-- if [[ $TRAVIS_PYTHON_VERSION == "3.3" ]]; then pip install "sphinx<1.5.0"; fi
-- if [[ $TRAVIS_PYTHON_VERSION == "3.4" ]]; then pip install "matplotlib>=1.5.0"; fi
+- if [[ $TRAVIS_PYTHON_VERSION == "3.4" ]]; then pip install "matplotlib>=1.5.0,!=2.1.0"; fi
- if [[ $TRAVIS_PYTHON_VERSION == "3.4" ]]; then pip install "sphinx>=1.5.0"; fi
-- if [[ $TRAVIS_PYTHON_VERSION == "3.5" ]]; then pip install "matplotlib>=1.5.0"; fi
+- if [[ $TRAVIS_PYTHON_VERSION == "3.5" ]]; then pip install "matplotlib>=1.5.0,!=2.1.0"; fi
- if [[ $TRAVIS_PYTHON_VERSION == "3.5" ]]; then pip install "sphinx>=1.5.0"; fi
-- if [[ $TRAVIS_PYTHON_VERSION == "3.6" ]]; then pip install "matplotlib>=1.5.0"; fi
+- if [[ $TRAVIS_PYTHON_VERSION == "3.6" ]]; then pip install "matplotlib>=1.5.0,!=2.1.0"; fi
- if [[ $TRAVIS_PYTHON_VERSION == "3.6" ]]; then pip install "sphinx>=1.5.0"; fi
- pip install -r requirements.txt
- pip install -e ".[pykdtree,quicklook]"
- pip install coveralls
script:
-# Once doctests pass this should be uncommented to replace the other coverage run line
- coverage run --source=pyresample setup.py test && cd docs && mkdir doctest && sphinx-build -E -n -b doctest ./source ./doctest && cd ..
after_success: coveralls
notifications:
diff --git a/appveyor.yml b/appveyor.yml
index ef6de64..51f60d1 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -28,13 +28,13 @@ environment:
PYTHON_ARCH: "64"
MINICONDA_VERSION: "3"
- - PYTHON: "C:\\Python35_32"
- PYTHON_VERSION: "3.5"
+ - PYTHON: "C:\\Python36_32"
+ PYTHON_VERSION: "3.6"
PYTHON_ARCH: "32"
MINICONDA_VERSION: "3"
- - PYTHON: "C:\\Python35_64"
- PYTHON_VERSION: "3.5"
+ - PYTHON: "C:\\Python36_64"
+ PYTHON_VERSION: "3.6"
PYTHON_ARCH: "64"
MINICONDA_VERSION: "3"
@@ -88,8 +88,3 @@ artifacts:
#on_success:
# - TODO: upload the content of dist/*.whl to a public wheelhouse
#
-
-notifications:
- - provider: Slack
- incoming_webhook:
- secure: 7bEOYCIxHE5PkCF156zjxbJIeKkv7UpZulyn+Di09jqDlpvZjR0Qj3a1LK9AOjgwWgaQYbeI4BEYEdeq7pPxDs9snK8qvvbFDbuLAzg+HEQ=
diff --git a/changelog.rst b/changelog.rst
index 820a32d..090290a 100644
--- a/changelog.rst
+++ b/changelog.rst
@@ -2,6 +2,103 @@ Changelog
=========
+v1.7.0 (2017-10-13)
+-------------------
+- update changelog. [Martin Raspaud]
+- Bump version: 1.6.1 → 1.7.0. [Martin Raspaud]
+- Merge pull request #82 from pytroll/fix-resample-bilinear. [David
+ Hoese]
+
+ Fix output shape of resample_bilinear()
+- Reshape output to have correct shape for the output area and num of
+ chans. [Panu Lahtinen]
+- Update tests to check proper output shape for resample_bilinear()
+ [Panu Lahtinen]
+- Merge pull request #79 from pytroll/fix-bil-documentation. [David
+ Hoese]
+
+ Fix example data for BIL, clarify text and add missing output_shape p…
+- Merge branch 'fix-bil-documentation' of
+ https://github.com/mraspaud/pyresample into fix-bil-documentation.
+ [Panu Lahtinen]
+- Fix example data for BIL, clarify text and add missing output_shape
+ param. [Panu Lahtinen]
+- Fix example data for BIL, clarify text and add missing output_shape
+ param. [Panu Lahtinen]
+- Merge pull request #75 from pytroll/fix-bil-mask-deprecation. [David
+ Hoese]
+
+ Fix bil mask deprecation
+- Merge branch 'develop' into fix-bil-mask-deprecation. [David Hoese]
+- Merge pull request #81 from pytroll/fix-reduce-bil-memory-use. [David
+ Hoese]
+
+ Reduce the memory use for ImageContainerBilinear tests
+- Reduce area size for BIL, reduce neighbours and adjust expected
+ results. [Panu Lahtinen]
+- Add proj4_dict_to_str utility function (#78) [David Hoese]
+
+ * Add proj4_dict_to_str utility function
+
+ Includes fixes for dynamic area definitions proj_id and
+ small performance improvement for projection coordinate generation
+
+ * Use more descriptive variable names
+
+ * Fix proj4 dict conversion test
+
+ * Exclude buggy version of matplotlib in travis tests
+
+ * Change appveyor python 3.5 environments to python 3.6
+
+ Also removes slack notification webhook which is no longer the
+ recommended way to post to slack from appveyor.
+
+ * Fix proj4 dict to string against recent changes to str to dict funcs
+
+- Utils edits for retreiving projection semi-major / semi-minor axes
+ (#77) [goodsonr]
+
+ proj4 strings converted to dictionary now consistent with other code (no longer has leading '+')
+ new logic for reporting projection semi-major / semi-minor axes ('a', 'b') based on information in proj4
+
+- Merge pull request #71 from pytroll/feature-bilinear-image. [David
+ Hoese]
+
+ Add image container for bilinear interpolation
+- Fix test result assertation. [Panu Lahtinen]
+- Add tests for ImageContainerBilinear, rewrap long lines. [Panu
+ Lahtinen]
+- Fix docstrings. [Panu Lahtinen]
+- Mention also ImageContainerBilinear. [Panu Lahtinen]
+- Handle 3D input data with bilinear interpolation. [Panu Lahtinen]
+- Add ImageContainerBilinear, autopep8. [Panu Lahtinen]
+- Merge pull request #74 from pytroll/fix-close-area-file. [David Hoese]
+
+ Use context manager to open area definition files
+- Use context manager to open files, PEP8. [Panu Lahtinen]
+- Merge pull request #76 from pytroll/feature-xarray. [Martin Raspaud]
+
+ Support resampling of xarray.DataArrays
+- Move docstring to init for consistency. [Martin Raspaud]
+- Merge develop into feature_xarray. [Martin Raspaud]
+- Support get_lonlats_dask in StackedAreaDefinitions. [Martin Raspaud]
+- Add get_lonlats_dask for SwathDefinitions. [Martin Raspaud]
+- Fix resampling of multidimensional xarrays. [Martin Raspaud]
+- Support xarray and use dask for simple cases. [Martin Raspaud]
+- WIP: Resampler for xarrays using dask. [Martin Raspaud]
+- Fix formatting. [Martin Raspaud]
+- Optimize memory consumption. [Martin Raspaud]
+- Clean up doc formatting. [Martin Raspaud]
+- Add dask.Array returning get_lonlats and get_proj_coords. [Martin
+ Raspaud]
+- Remove Python 3.3 from travis tests, it's not supported anymore. [Panu
+ Lahtinen]
+- Supress UserWarning about possible extra neighbours within search
+ radius. [Panu Lahtinen]
+- Handle masked arrays properly for new Numpy versions. [Panu Lahtinen]
+
+
v1.6.1 (2017-09-18)
-------------------
- update changelog. [Martin Raspaud]
diff --git a/docs/source/swath.rst b/docs/source/swath.rst
index b370996..0e06930 100644
--- a/docs/source/swath.rst
+++ b/docs/source/swath.rst
@@ -8,7 +8,7 @@ Resampling can be done using nearest neighbour method, Guassian weighting, weigh
pyresample.image
----------------
-The ImageContainerNearest class can be used for nearest neighbour resampling of swaths as well as grids.
+The ImageContainerNearest and ImageContanerBilinear classes can be used for resampling of swaths as well as grids. Below is an example using nearest neighbour resampling.
.. doctest::
@@ -29,7 +29,7 @@ The ImageContainerNearest class can be used for nearest neighbour resampling of
>>> area_con = swath_con.resample(area_def)
>>> result = area_con.image_data
-For other resampling types or splitting the process in two steps use the functions in **pyresample.kd_tree** described below.
+For other resampling types or splitting the process in two steps use e.g. the functions in **pyresample.kd_tree** described below.
pyresample.kd_tree
------------------
@@ -282,9 +282,9 @@ Function for resampling using bilinear interpolation for irregular source grids.
... 800, 800,
... [-1370912.72, -909968.64,
... 1029087.28, 1490031.36])
- >>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
- >>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
- >>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
+ >>> data = np.fromfunction(lambda y, x: y*x, (500, 100))
+ >>> lons = np.fromfunction(lambda y, x: 3 + x * 0.1, (500, 100))
+ >>> lats = np.fromfunction(lambda y, x: 75 - y * 0.1, (500, 100))
>>> source_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> result = bilinear.resample_bilinear(data, source_def, target_def,
... radius=50e3, neighbours=32,
@@ -304,7 +304,8 @@ Keyword arguments which are passed to **kd_tree**:
* **radius**: radius around each target pixel in meters to search for
neighbours in the source data
* **neighbours**: number of closest locations to consider when
- selecting the four data points around the target pixel
+ selecting the four data points around the target location. Note that this
+ value needs to be large enough to ensure "surrounding" the target!
* **nprocs**: number of processors to use for finding the closest pixels
* **fill_value**: fill invalid pixel with this value. If
**fill_value=None** is used, masked arrays will be returned
@@ -332,7 +333,9 @@ significantly. This is also done internally by the
**resample_bilinear** function, but separating these steps makes it
possible to cache the coefficients if the same transformation is done
over and over again. This is very typical in operational
-geostationary satellite image processing.
+geostationary satellite image processing. Note that the output shape is now
+defined so that the result is reshaped to correct shape. This reshaping
+is done internally in **resample_bilinear**.
.. doctest::
@@ -353,7 +356,8 @@ geostationary satellite image processing.
>>> t_params, s_params, input_idxs, idx_ref = \
... bilinear.get_bil_info(source_def, target_def, radius=50e3, nprocs=1)
>>> res = bilinear.get_sample_from_bil_info(data.ravel(), t_params, s_params,
- ... input_idxs, idx_ref)
+ ... input_idxs, idx_ref,
+ ... output_shape=target_def.shape)
pyresample.ewa
diff --git a/pyresample/bilinear/__init__.py b/pyresample/bilinear/__init__.py
index 4f721f1..89f7655 100644
--- a/pyresample/bilinear/__init__.py
+++ b/pyresample/bilinear/__init__.py
@@ -30,6 +30,7 @@ http://www.ahinson.com/algorithms_general/Sections/InterpolationRegression/Inter
import numpy as np
from pyproj import Proj
+import warnings
from pyresample import kd_tree
@@ -97,6 +98,12 @@ def resample_bilinear(data, source_geo_def, target_area_def, radius=50e3,
else:
result[np.isnan(result)] = fill_value
+ # Reshape to target area shape
+ shp = target_area_def.shape
+ result = result.reshape((shp[0], shp[1], data.shape[1]))
+ # Remove extra dimensions
+ result = np.squeeze(result)
+
return result
@@ -210,11 +217,13 @@ def get_bil_info(source_geo_def, target_area_def, radius=50e3, neighbours=32,
# source_geo_def = SwathDefinition(lons, lats)
# Calculate neighbour information
- (input_idxs, output_idxs, idx_ref, dists) = \
- kd_tree.get_neighbour_info(source_geo_def, target_area_def,
- radius, neighbours=neighbours,
- nprocs=nprocs, reduce_data=reduce_data,
- segments=segments, epsilon=epsilon)
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore")
+ (input_idxs, output_idxs, idx_ref, dists) = \
+ kd_tree.get_neighbour_info(source_geo_def, target_area_def,
+ radius, neighbours=neighbours,
+ nprocs=nprocs, reduce_data=reduce_data,
+ segments=segments, epsilon=epsilon)
del output_idxs, dists
@@ -389,8 +398,14 @@ def _mask_coordinates(lons, lats):
lats = lats.ravel()
idxs = ((lons < -180.) | (lons > 180.) |
(lats < -90.) | (lats > 90.))
- lons[idxs] = np.nan
- lats[idxs] = np.nan
+ if hasattr(lons, 'mask'):
+ lons = np.ma.masked_where(idxs | lons.mask, lons)
+ else:
+ lons[idxs] = np.nan
+ if hasattr(lats, 'mask'):
+ lats = np.ma.masked_where(idxs | lats.mask, lats)
+ else:
+ lats[idxs] = np.nan
return lons, lats
diff --git a/pyresample/geometry.py b/pyresample/geometry.py
index 31d84ee..988c310 100644
--- a/pyresample/geometry.py
+++ b/pyresample/geometry.py
@@ -465,6 +465,14 @@ class SwathDefinition(CoordinateDefinition):
elif lons.ndim > 2:
raise ValueError('Only 1 and 2 dimensional swaths are allowed')
+ def get_lonlats_dask(self, blocksize=1000, dtype=None):
+ """Get the lon lats as a single dask array."""
+ import dask.array as da
+ lons, lats = self.get_lonlats()
+ lons = da.from_array(lons, chunks=blocksize * lons.ndim)
+ lats = da.from_array(lats, chunks=blocksize * lons.ndim)
+ return da.stack((lons, lats), axis=-1)
+
def _compute_omerc_parameters(self, ellipsoid):
"""Compute the oblique mercator projection bouding box parameters."""
lines, cols = self.lons.shape
@@ -622,7 +630,7 @@ class DynamicAreaDefinition(object):
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,
+ return AreaDefinition(self.area_id, self.description, '',
self.proj_dict, self.x_size, self.y_size,
self.area_extent)
@@ -751,6 +759,10 @@ class AreaDefinition(BaseDefinition):
self.dtype = dtype
+ @property
+ def proj_str(self):
+ return utils.proj4_dict_to_str(self.proj_dict, sort=True)
+
def __str__(self):
# We need a sorted dictionary for a unique hash of str(self)
proj_dict = self.proj_dict
@@ -759,7 +771,7 @@ class AreaDefinition(BaseDefinition):
for k in sorted(proj_dict.keys())]) +
'}')
- if self.proj_id is None:
+ if not self.proj_id:
third_line = ""
else:
third_line = "Projection ID: {0}\n".format(self.proj_id)
@@ -951,6 +963,34 @@ class AreaDefinition(BaseDefinition):
return self.get_lonlats(nprocs=None, data_slice=(row, col))
+ def get_proj_coords_dask(self, blocksize, dtype=None):
+ from dask.base import tokenize
+ from dask.array import Array
+ if dtype is None:
+ dtype = self.dtype
+
+ vchunks = range(0, self.y_size, blocksize)
+ hchunks = range(0, self.x_size, blocksize)
+
+ token = tokenize(blocksize)
+ name = 'get_proj_coords-' + token
+
+ dskx = {(name, i, j, 0): (self.get_proj_coords_array,
+ (slice(vcs, min(vcs + blocksize, self.y_size)),
+ slice(hcs, min(hcs + blocksize, self.x_size))),
+ False, dtype)
+ for i, vcs in enumerate(vchunks)
+ for j, hcs in enumerate(hchunks)
+ }
+
+ res = Array(dskx, name, shape=list(self.shape) + [2],
+ chunks=(blocksize, blocksize, 2),
+ dtype=dtype)
+ return res
+
+ def get_proj_coords_array(self, data_slice=None, cache=False, dtype=None):
+ return np.dstack(self.get_proj_coords(data_slice, cache, dtype))
+
def get_proj_coords(self, data_slice=None, cache=False, dtype=None):
"""Get projection coordinates of grid
@@ -959,7 +999,7 @@ class AreaDefinition(BaseDefinition):
data_slice : slice object, optional
Calculate only coordinates for specified slice
cache : bool, optional
- Store result the result. Requires data_slice to be None
+ Store the result. Requires data_slice to be None
Returns
-------
@@ -1041,17 +1081,9 @@ class AreaDefinition(BaseDefinition):
cols = 1
# Calculate coordinates
- target_x = np.fromfunction(lambda i, j: (j + col_start) *
- self.pixel_size_x +
- self.pixel_upper_left[0],
- (rows,
- cols), dtype=dtype)
-
- target_y = np.fromfunction(lambda i, j:
- self.pixel_upper_left[1] -
- (i + row_start) * self.pixel_size_y,
- (rows,
- cols), dtype=dtype)
+ target_x = np.arange(col_start, col_start + cols, dtype=dtype) * self.pixel_size_x + self.pixel_upper_left[0]
+ target_y = np.arange(row_start, row_start + rows, dtype=dtype) * -self.pixel_size_y + self.pixel_upper_left[1]
+ target_x, target_y = np.meshgrid(target_x, target_y)
if is_single_value:
# Return single values
@@ -1079,12 +1111,14 @@ class AreaDefinition(BaseDefinition):
@property
def proj_x_coords(self):
- warnings.warn("Deprecated, use 'projection_x_coords' instead", DeprecationWarning)
+ warnings.warn(
+ "Deprecated, use 'projection_x_coords' instead", DeprecationWarning)
return self.projection_x_coords
@property
def proj_y_coords(self):
- warnings.warn("Deprecated, use 'projection_y_coords' instead", DeprecationWarning)
+ warnings.warn(
+ "Deprecated, use 'projection_y_coords' instead", DeprecationWarning)
return self.projection_y_coords
@property
@@ -1104,6 +1138,39 @@ class AreaDefinition(BaseDefinition):
Coordinate(corner_lons[2], corner_lats[2]),
Coordinate(corner_lons[3], corner_lats[3])]
+ def get_lonlats_dask(self, blocksize=1000, dtype=None):
+ from dask.base import tokenize
+ from dask.array import Array
+ import pyproj
+
+ dtype = dtype or self.dtype
+ proj_coords = self.get_proj_coords_dask(blocksize, dtype)
+ target_x, target_y = proj_coords[:, :, 0], proj_coords[:, :, 1]
+
+ target_proj = pyproj.Proj(**self.proj_dict)
+
+ def invproj(data1, data2):
+ return np.dstack(target_proj(data1.compute(), data2.compute(), inverse=True))
+ token = tokenize(str(self), blocksize, dtype)
+ name = 'get_lonlats-' + token
+
+ vchunks = range(0, self.y_size, blocksize)
+ hchunks = range(0, self.x_size, blocksize)
+
+ dsk = {(name, i, j, 0): (invproj,
+ target_x[slice(vcs, min(vcs + blocksize, self.y_size)),
+ slice(hcs, min(hcs + blocksize, self.x_size))],
+ target_y[slice(vcs, min(vcs + blocksize, self.y_size)),
+ slice(hcs, min(hcs + blocksize, self.x_size))])
+ for i, vcs in enumerate(vchunks)
+ for j, hcs in enumerate(hchunks)
+ }
+
+ res = Array(dsk, name, shape=list(self.shape) + [2],
+ chunks=(blocksize, blocksize, 2),
+ dtype=dtype)
+ return res
+
def get_lonlats(self, nprocs=None, data_slice=None, cache=False, dtype=None):
"""Returns lon and lat arrays of area.
@@ -1293,6 +1360,20 @@ class StackedAreaDefinition(BaseDefinition):
return self.lons, self.lats
+ def get_lonlats_dask(self, blocksize=1000, dtype=None):
+ """"Return lon and lat dask arrays of the area."""
+ import dask.array as da
+ llonslats = []
+ for definition in self.defs:
+ lonslats = definition.get_lonlats_dask(blocksize=blocksize,
+ dtype=dtype)
+
+ llonslats.append(lonslats)
+
+ self.lonlats = da.concatenate(llonslats, axis=0)
+
+ return self.lonlats
+
def squeeze(self):
"""Generate a single AreaDefinition if possible."""
if len(self.defs) == 1:
diff --git a/pyresample/image.py b/pyresample/image.py
index 65d7462..5254ac1 100644
--- a/pyresample/image.py
+++ b/pyresample/image.py
@@ -21,32 +21,32 @@ from __future__ import absolute_import
import numpy as np
-from pyresample import geometry, grid, kd_tree
+from pyresample import geometry, grid, kd_tree, bilinear
class ImageContainer(object):
- """Holds image with geometry definition.
+ """Holds image with geometry definition.
Allows indexing with linesample arrays.
Parameters
----------
image_data : numpy array
Image data
- geo_def : object
+ geo_def : object
Geometry definition
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
- nprocs : int, optional
+ nprocs : int, optional
Number of processor cores to be used
Attributes
----------
image_data : numpy array
Image data
- geo_def : object
+ geo_def : object
Geometry definition
fill_value : int or None
Resample result fill value
@@ -99,8 +99,8 @@ class ImageContainer(object):
----------
row_indices : numpy array
Row indices. Dimensions must match col_indices
- col_indices : numpy array
- Col indices. Dimensions must match row_indices
+ col_indices : numpy array
+ Col indices. Dimensions must match row_indices
Returns
-------
@@ -133,13 +133,13 @@ class ImageContainerQuick(ImageContainer):
----------
image_data : numpy array
Image data
- geo_def : object
+ geo_def : object
Area definition as AreaDefinition object
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
- nprocs : int, optional
+ nprocs : int, optional
Number of processor cores to be used for geometry operations
segments : int or None
Number of segments to use when resampling.
@@ -149,16 +149,16 @@ class ImageContainerQuick(ImageContainer):
----------
image_data : numpy array
Image data
- geo_def : object
+ geo_def : object
Area definition as AreaDefinition object
fill_value : int or None
Resample result fill 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
nprocs : int
Number of processor cores to be used
segments : int or None
- Number of segments to use when resampling
+ Number of segments to use when resampling
"""
def __init__(self, image_data, geo_def, fill_value=0, nprocs=1,
@@ -172,7 +172,7 @@ class ImageContainerQuick(ImageContainer):
self.segments = segments
def resample(self, target_area_def):
- """Resamples image to area definition using nearest neighbour
+ """Resamples image to area definition using nearest neighbour
approach in projection coordinates.
Parameters
@@ -183,7 +183,7 @@ class ImageContainerQuick(ImageContainer):
Returns
-------
image_container : object
- ImageContainerQuick object of resampled area
+ ImageContainerQuick object of resampled area
"""
resampled_image = grid.get_resampled_image(target_area_def,
@@ -200,28 +200,28 @@ class ImageContainerQuick(ImageContainer):
class ImageContainerNearest(ImageContainer):
- """Holds image with geometry definition.
- Allows nearest neighbour resampling to new geometry definition.
+ """Holds image with geometry definition.
+ Allows nearest neighbour to new geometry definition.
Parameters
----------
image_data : numpy array
Image data
- geo_def : object
+ geo_def : object
Geometry definition
- radius_of_influence : float
- Cut off distance in meters
+ 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
+ If fill_value is None a masked array is returned
with undetermined pixels masked
reduce_data : bool, optional
Perform coarse data reduction before resampling in order
to reduce execution time
- nprocs : int, optional
+ nprocs : int, optional
Number of processor cores to be used for geometry operations
segments : int or None
Number of segments to use when resampling.
@@ -230,12 +230,12 @@ class ImageContainerNearest(ImageContainer):
Attributes
----------
- image_data : numpy array
+ image_data : numpy array
Image data
- geo_def : object
+ geo_def : object
Geometry definition
- radius_of_influence : float
- Cut off distance in meters
+ radius_of_influence : float
+ Cut off distance in meters
epsilon : float
Allowed uncertainty in meters
fill_value : int or None
@@ -245,7 +245,7 @@ class ImageContainerNearest(ImageContainer):
nprocs : int
Number of processor cores to be used
segments : int or None
- Number of segments to use when resampling
+ Number of segments to use when resampling
"""
def __init__(self, image_data, geo_def, radius_of_influence, epsilon=0,
@@ -259,18 +259,18 @@ class ImageContainerNearest(ImageContainer):
self.segments = segments
def resample(self, target_geo_def):
- """Resamples image to area definition using nearest neighbour
+ """Resamples image to area definition using nearest neighbour
approach
Parameters
----------
- target_geo_def : object
- Target geometry definition
+ target_geo_def : object
+ Target geometry definition
Returns
-------
image_container : object
- ImageContainerNearest object of resampled geometry
+ ImageContainerNearest object of resampled geometry
"""
if self.image_data.ndim > 2 and self.ndim > 1:
@@ -297,3 +297,121 @@ class ImageContainerNearest(ImageContainer):
reduce_data=self.reduce_data,
nprocs=self.nprocs,
segments=self.segments)
+
+
+class ImageContainerBilinear(ImageContainer):
+
+ """Holds image with geometry definition.
+ Allows bilinear to new geometry definition.
+
+ Parameters
+ ----------
+ image_data : numpy array
+ Image data
+ geo_def : object
+ Geometry definition
+ 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
+ reduce_data : bool, optional
+ Perform coarse data reduction before resampling in order
+ to reduce execution time
+ nprocs : int, optional
+ Number of processor cores to be used for geometry operations
+ segments : int or None
+ Number of segments to use when resampling.
+ If set to None an estimate will be calculated
+
+ Attributes
+ ----------
+
+ image_data : numpy array
+ Image data
+ geo_def : object
+ Geometry definition
+ radius_of_influence : float
+ Cut off distance in meters
+ epsilon : float
+ Allowed uncertainty in meters
+ fill_value : int or None
+ Resample result fill value
+ reduce_data : bool
+ Perform coarse data reduction before resampling
+ nprocs : int
+ Number of processor cores to be used
+ segments : int or None
+ Number of segments to use when resampling
+ """
+
+ def __init__(self, image_data, geo_def, radius_of_influence, epsilon=0,
+ fill_value=0, reduce_data=True, nprocs=1, segments=None,
+ neighbours=32):
+ super(ImageContainerBilinear, self).__init__(image_data, geo_def,
+ fill_value=fill_value,
+ nprocs=nprocs)
+ self.radius_of_influence = radius_of_influence
+ self.epsilon = epsilon
+ self.reduce_data = reduce_data
+ self.segments = segments
+ self.neighbours = neighbours
+
+ def resample(self, target_geo_def):
+ """Resamples image to area definition using bilinear approach
+
+ Parameters
+ ----------
+ target_geo_def : object
+ Target geometry definition
+
+ Returns
+ -------
+ image_container : object
+ ImageContainerBilinear object of resampled geometry
+ """
+
+ if self.image_data.ndim > 2 and self.ndim > 1:
+ image_data = self.image_data.reshape(self.image_data.shape[0] *
+ self.image_data.shape[1],
+ self.image_data.shape[2])
+ else:
+ image_data = self.image_data.ravel()
+
+ try:
+ mask = image_data.mask.copy()
+ image_data = image_data.data.copy()
+ image_data[mask] = np.nan
+ except AttributeError:
+ pass
+
+ resampled_image = \
+ bilinear.resample_bilinear(image_data,
+ self.geo_def,
+ target_geo_def,
+ radius=self.radius_of_influence,
+ neighbours=self.neighbours,
+ epsilon=self.epsilon,
+ fill_value=self.fill_value,
+ nprocs=self.nprocs,
+ reduce_data=self.reduce_data,
+ segments=self.segments)
+ try:
+ resampled_image = resampled_image.reshape(target_geo_def.shape)
+ except ValueError:
+ # The input data was 3D
+ shp = target_geo_def.shape
+ new_shp = [shp[0], shp[1], image_data.shape[-1]]
+ resampled_image = resampled_image.reshape(new_shp)
+
+ return ImageContainerBilinear(resampled_image, target_geo_def,
+ self.radius_of_influence,
+ epsilon=self.epsilon,
+ fill_value=self.fill_value,
+ reduce_data=self.reduce_data,
+ nprocs=self.nprocs,
+ segments=self.segments)
diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py
index e8cdcd2..a838461 100644
--- a/pyresample/kd_tree.py
+++ b/pyresample/kd_tree.py
@@ -24,11 +24,20 @@ from __future__ import absolute_import
import sys
import types
import warnings
+from logging import getLogger
import numpy as np
from pyresample import _spatial_mp, data_reduce, geometry
+logger = getLogger(__name__)
+
+try:
+ from xarray import DataArray
+ import dask.array as da
+except ImportError:
+ logger.info("XArray or dask unavailable, some functionality missing.")
+
if sys.version < '3':
range = xrange
else:
@@ -877,6 +886,274 @@ def get_sample_from_neighbour_info(resample_type, output_shape, data,
return result
+class XArrayResamplerNN(object):
+
+ def __init__(self, source_geo_def, target_geo_def, radius_of_influence,
+ neighbours=8, epsilon=0, reduce_data=True,
+ nprocs=1, segments=None):
+ """
+ Parameters
+ ----------
+ source_geo_def : object
+ Geometry definition of source
+ target_geo_def : object
+ Geometry definition of target
+ radius_of_influence : float
+ Cut off distance in meters
+ 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
+ reduce_data : bool, optional
+ Perform initial coarse reduction of source dataset in order
+ to reduce execution time
+ nprocs : int, optional
+ Number of processor cores to be used
+ segments : int or None
+ Number of segments to use when resampling.
+ If set to None an estimate will be calculated
+ """
+
+ self.valid_input_index = None
+ self.valid_output_index = None
+ self.index_array = None
+ self.distance_array = None
+ self.neighbours = neighbours
+ self.epsilon = epsilon
+ self.reduce_data = reduce_data
+ self.nprocs = nprocs
+ self.segments = segments
+ self.source_geo_def = source_geo_def
+ self.target_geo_def = target_geo_def
+ self.radius_of_influence = radius_of_influence
+
+ def transform_lonlats(self, 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))
+
+ return da.stack((x_coords, y_coords, z_coords), axis=-1)
+
+ def _create_resample_kdtree(self, source_lons, source_lats, valid_input_index):
+ """Set up kd tree on input"""
+
+ """
+ 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)
+ input_coords = input_coords[valid_input_index]
+ """
+
+ vii = valid_input_index.compute().ravel()
+ source_lons_valid = source_lons.ravel()[vii]
+ source_lats_valid = source_lats.ravel()[vii]
+
+ input_coords = self.transform_lonlats(source_lons_valid,
+ source_lats_valid)
+
+ if input_coords.size == 0:
+ raise EmptyResult('No valid data points in input data')
+
+ # Build kd-tree on input
+
+ if kd_tree_name == 'pykdtree':
+ resample_kdtree = KDTree(input_coords.compute())
+ else:
+ resample_kdtree = sp.cKDTree(input_coords.compute())
+
+ return resample_kdtree
+
+ def _query_resample_kdtree(self, resample_kdtree, target_lons,
+ target_lats, valid_output_index,
+ reduce_data=True):
+ """Query kd-tree on slice of target coordinates"""
+ from dask.base import tokenize
+ from dask.array import Array
+
+ def query(target_lons, target_lats, valid_output_index, c_slice):
+ voi = valid_output_index[c_slice].compute()
+ shape = voi.shape
+ voir = voi.ravel()
+ target_lons_valid = target_lons[c_slice].ravel()[voir]
+ target_lats_valid = target_lats[c_slice].ravel()[voir]
+
+ coords = self.transform_lonlats(target_lons_valid,
+ target_lats_valid)
+ distance_array, index_array = np.stack(
+ resample_kdtree.query(coords.compute(),
+ k=self.neighbours,
+ eps=self.epsilon,
+ distance_upper_bound=self.radius_of_influence))
+
+ res_ia = np.full(shape, fill_value=np.nan, dtype=np.float)
+ res_da = np.full(shape, fill_value=np.nan, dtype=np.float)
+ res_ia[voi] = index_array
+ res_da[voi] = distance_array
+ return np.stack([res_ia, res_da], axis=-1)
+
+ token = tokenize(1000)
+ name = 'query-' + token
+
+ dsk = {}
+ vstart = 0
+
+ for i, vck in enumerate(valid_output_index.chunks[0]):
+ hstart = 0
+ for j, hck in enumerate(valid_output_index.chunks[1]):
+ c_slice = (slice(vstart, vstart + vck),
+ slice(hstart, hstart + hck))
+ dsk[(name, i, j, 0)] = (query, target_lons,
+ target_lats, valid_output_index,
+ c_slice)
+ hstart += hck
+ vstart += vck
+
+ res = Array(dsk, name,
+ shape=list(valid_output_index.shape) + [2],
+ chunks=list(valid_output_index.chunks) + [2],
+ dtype=target_lons.dtype)
+
+ index_array = res[:, :, 0].astype(np.uint)
+ distance_array = res[:, :, 1]
+ return index_array, distance_array
+
+ def get_neighbour_info(self):
+ """Returns neighbour info
+
+ Returns
+ -------
+ (valid_input_index, valid_output_index,
+ index_array, distance_array) : tuple of numpy arrays
+ Neighbour resampling info
+ """
+
+ if self.source_geo_def.size < self.neighbours:
+ warnings.warn('Searching for %s neighbours in %s data points' %
+ (self.neighbours, self.source_geo_def.size))
+
+ source_lonlats = self.source_geo_def.get_lonlats_dask()
+ source_lons = source_lonlats[:, :, 0]
+ source_lats = source_lonlats[:, :, 1]
+ valid_input_index = ((source_lons >= -180) & (source_lons <= 180) &
+ (source_lats <= 90) & (source_lats >= -90))
+
+ # Create kd-tree
+ try:
+ resample_kdtree = self._create_resample_kdtree(source_lons,
+ source_lats,
+ valid_input_index)
+
+ except EmptyResult:
+ # Handle if all input data is reduced away
+ valid_output_index, index_array, distance_array = \
+ _create_empty_info(self.source_geo_def,
+ self.target_geo_def, self.neighbours)
+ return (valid_input_index, valid_output_index, index_array,
+ distance_array)
+
+ target_lonlats = self.target_geo_def.get_lonlats_dask()
+ target_lons = target_lonlats[:, :, 0]
+ target_lats = target_lonlats[:, :, 1]
+ valid_output_index = ((target_lons >= -180) & (target_lons <= 180) &
+ (target_lats <= 90) & (target_lats >= -90))
+
+ index_array, distance_array = self._query_resample_kdtree(resample_kdtree,
+ target_lons,
+ target_lats,
+ valid_output_index)
+
+ self.valid_input_index = valid_input_index
+ self.valid_output_index = valid_output_index
+ self.index_array = index_array
+ self.distance_array = distance_array
+
+ return valid_input_index, valid_output_index, index_array, distance_array
+
+ def get_sample_from_neighbour_info(self, data):
+
+ # flatten x and y in the source array
+
+ output_shape = []
+ chunks = []
+ source_dims = data.dims
+ for dim in source_dims:
+ if dim == 'y':
+ output_shape += [self.target_geo_def.y_size]
+ chunks += [1000]
+ elif dim == 'x':
+ output_shape += [self.target_geo_def.x_size]
+ chunks += [1000]
+ else:
+ output_shape += [data[dim].size]
+ chunks += [10]
+
+ new_dims = []
+ xy_dims = []
+ source_shape = [1, 1]
+ chunks = [1, 1]
+ for i, dim in enumerate(data.dims):
+ if dim not in ['x', 'y']:
+ new_dims.append(dim)
+ source_shape[1] *= data.shape[i]
+ chunks[1] *= 10
+ else:
+ xy_dims.append(dim)
+ source_shape[0] *= data.shape[i]
+ chunks[0] *= 1000
+
+ new_dims = xy_dims + new_dims
+
+ target_shape = [np.prod(self.target_geo_def.shape), source_shape[1]]
+ source_data = data.transpose(*new_dims).data.reshape(source_shape)
+
+ input_size = self.valid_input_index.sum()
+ index_mask = (self.index_array == input_size)
+ new_index_array = da.where(
+ index_mask, 0, self.index_array).ravel().astype(int).compute()
+ valid_targets = self.valid_output_index.ravel()
+
+ target_lines = []
+
+ for line in range(target_shape[1]):
+ #target_data_line = target_data[:, line]
+ new_data = source_data[:, line][self.valid_input_index.ravel()]
+ # could this be a bug in dask ? we have to compute to avoid errors
+ result = new_data.compute()[new_index_array]
+ result[index_mask.ravel()] = np.nan
+ #target_data_line = da.full(target_shape[0], np.nan, chunks=1000000)
+ target_data_line = np.full(target_shape[0], np.nan)
+ target_data_line[valid_targets] = result
+ target_lines.append(target_data_line[:, np.newaxis])
+
+ target_data = np.hstack(target_lines)
+
+ new_shape = []
+ for dim in new_dims:
+ if dim == 'x':
+ new_shape.append(self.target_geo_def.x_size)
+ elif dim == 'y':
+ new_shape.append(self.target_geo_def.y_size)
+ else:
+ new_shape.append(data[dim].size)
+
+ output_arr = DataArray(da.from_array(target_data.reshape(new_shape), chunks=[1000] * len(new_shape)),
+ dims=new_dims)
+ for dim in source_dims:
+ if dim == 'x':
+ output_arr['x'] = self.target_geo_def.proj_x_coords
+ elif dim == 'y':
+ output_arr['y'] = self.target_geo_def.proj_y_coords
+ else:
+ output_arr[dim] = data[dim]
+
+ return output_arr.transpose(*source_dims)
+
+
def _get_fill_mask_value(data_dtype):
"""Returns the maximum value of dtype"""
diff --git a/pyresample/test/test_bilinear.py b/pyresample/test/test_bilinear.py
index 7bb3362..7a0eb1c 100644
--- a/pyresample/test/test_bilinear.py
+++ b/pyresample/test/test_bilinear.py
@@ -207,7 +207,7 @@ class Test(unittest.TestCase):
res = bil.resample_bilinear(self.data1,
self.swath_def,
self.target_def)
- self.assertEqual(res.size, self.target_def.size)
+ self.assertEqual(res.shape, self.target_def.shape)
# There should be only one pixel with value 1, all others are 0
self.assertEqual(res.sum(), 1)
@@ -225,8 +225,8 @@ class Test(unittest.TestCase):
self.swath_def,
self.target_def)
shp = res.shape
- self.assertEqual(shp[0], self.target_def.size)
- self.assertEqual(shp[1], 2)
+ self.assertEqual(shp[0:2], self.target_def.shape)
+ self.assertEqual(shp[-1], 2)
def suite():
diff --git a/pyresample/test/test_image.py b/pyresample/test/test_image.py
index 8342583..b6bbd4f 100644
--- a/pyresample/test/test_image.py
+++ b/pyresample/test/test_image.py
@@ -18,7 +18,8 @@ def tmp(f):
class Test(unittest.TestCase):
- area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
+ area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
+ 'areaD',
{'a': '6378144.0',
'b': '6356759.0',
'lat_0': '50.00',
@@ -28,11 +29,12 @@ class Test(unittest.TestCase):
800,
800,
[-1370912.72,
- -909968.64000000001,
- 1029087.28,
- 1490031.3600000001])
+ -909968.64000000001,
+ 1029087.28,
+ 1490031.3600000001])
- msg_area = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
+ msg_area = geometry.AreaDefinition('msg_full',
+ 'Full globe MSG image 0 degrees',
'msg_full',
{'a': '6378169.0',
'b': '6356584.0',
@@ -42,12 +44,12 @@ class Test(unittest.TestCase):
3712,
3712,
[-5568742.4000000004,
- -5568742.4000000004,
- 5568742.4000000004,
- 5568742.4000000004]
- )
+ -5568742.4000000004,
+ 5568742.4000000004,
+ 5568742.4000000004])
- msg_area_resize = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
+ msg_area_resize = geometry.AreaDefinition('msg_full',
+ 'Full globe MSG image 0 degrees',
'msg_full',
{'a': '6378169.0',
'b': '6356584.0',
@@ -57,10 +59,9 @@ class Test(unittest.TestCase):
928,
928,
[-5568742.4000000004,
- -5568742.4000000004,
- 5568742.4000000004,
- 5568742.4000000004]
- )
+ -5568742.4000000004,
+ 5568742.4000000004,
+ 5568742.4000000004])
@tmp
def test_image(self):
@@ -103,7 +104,8 @@ class Test(unittest.TestCase):
area_con = msg_con.resample(self.area_def)
res = area_con.image_data
resampled_mask = res.mask.astype('int')
- expected = numpy.fromfile(os.path.join(os.path.dirname(__file__), 'test_files', 'mask_grid.dat'),
+ expected = numpy.fromfile(os.path.join(os.path.dirname(__file__),
+ 'test_files', 'mask_grid.dat'),
sep=' ').reshape((800, 800))
self.assertTrue(numpy.array_equal(
resampled_mask, expected), msg='Failed to resample masked array')
@@ -119,7 +121,9 @@ class Test(unittest.TestCase):
area_con = msg_con.resample(self.area_def)
res = area_con.image_data
resampled_mask = res.mask.astype('int')
- expected = numpy.fromfile(os.path.join(os.path.dirname(__file__), 'test_files', 'mask_grid.dat'),
+ expected = numpy.fromfile(os.path.join(os.path.dirname(__file__),
+ 'test_files',
+ 'mask_grid.dat'),
sep=' ').reshape((800, 800))
self.assertTrue(numpy.array_equal(
resampled_mask, expected), msg='Failed to resample masked array')
@@ -170,7 +174,8 @@ class Test(unittest.TestCase):
lambda y, x: y * x * 10 ** -6, (3712, 3712)) * 2
data = numpy.dstack((data1, data2))
msg_con = image.ImageContainer(data, self.msg_area)
- #area_con = msg_con.resample_area_nearest_neighbour(self.area_def, 50000)
+ # area_con = msg_con.resample_area_nearest_neighbour(self.area_def,
+ # 50000)
row_indices, col_indices = \
utils.generate_nearest_neighbour_linesample_arrays(self.msg_area,
self.area_def,
@@ -214,6 +219,47 @@ class Test(unittest.TestCase):
self.assertEqual(cross_sum, expected,
msg='ImageContainer swath segments resampling nearest failed')
+ def test_bilinear(self):
+ data = numpy.fromfunction(lambda y, x: y * x * 10 ** -6, (928, 928))
+ msg_con = image.ImageContainerBilinear(data, self.msg_area_resize,
+ 50000, segments=1,
+ neighbours=8)
+ area_con = msg_con.resample(self.area_def)
+ res = area_con.image_data
+ cross_sum = res.sum()
+ expected = 24690.127073654239
+ self.assertAlmostEqual(cross_sum, expected)
+
+ def test_bilinear_multi(self):
+ data1 = numpy.fromfunction(lambda y, x: y * x * 10 ** -6, (928, 928))
+ data2 = numpy.fromfunction(lambda y, x: y * x * 10 ** -6,
+ (928, 928)) * 2
+ data = numpy.dstack((data1, data2))
+ msg_con = image.ImageContainerBilinear(data, self.msg_area_resize,
+ 50000, segments=1,
+ neighbours=8)
+ area_con = msg_con.resample(self.area_def)
+ res = area_con.image_data
+ cross_sum1 = res[:, :, 0].sum()
+ expected1 = 24690.127073654239
+ self.assertAlmostEqual(cross_sum1, expected1)
+ cross_sum2 = res[:, :, 1].sum()
+ expected2 = 24690.127073654239 * 2
+ self.assertAlmostEqual(cross_sum2, expected2)
+
+ def test_bilinear_swath(self):
+ data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
+ lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
+ lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
+ swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
+ swath_con = image.ImageContainerBilinear(data, swath_def, 500000,
+ segments=1, neighbours=8)
+ area_con = swath_con.resample(self.area_def)
+ res = area_con.image_data
+ cross_sum = res.sum()
+ expected = 16762584.12441789
+ self.assertAlmostEqual(cross_sum, expected)
+
def suite():
"""The test suite.
diff --git a/pyresample/test/test_kd_tree.py b/pyresample/test/test_kd_tree.py
index eafc689..1ec401a 100644
--- a/pyresample/test/test_kd_tree.py
+++ b/pyresample/test/test_kd_tree.py
@@ -4,8 +4,9 @@ import os
import sys
import numpy
+
+from pyresample import data_reduce, geometry, kd_tree, utils
from pyresample.test.utils import catch_warnings
-from pyresample import kd_tree, utils, geometry, data_reduce
if sys.version_info < (2, 7):
import unittest2 as unittest
@@ -78,7 +79,7 @@ class Test(unittest.TestCase):
self.assertTrue(
len(w) > 0, 'Failed to create neighbour warning')
self.assertTrue((any('Searching' in str(_w.message) for _w in w)),
- 'Failed to create correct neighbour warning')
+ 'Failed to create correct neighbour warning')
expected_res = 2.20206560694
expected_stddev = 0.707115076173
@@ -101,7 +102,7 @@ class Test(unittest.TestCase):
self.assertTrue(
len(w) > 0, 'Failed to create neighbour warning')
self.assertTrue((any('Searching' in str(_w.message) for _w in w)),
- 'Failed to create correct neighbour warning')
+ 'Failed to create correct neighbour warning')
self.assertAlmostEqual(res[0], 2.32193149, 5,
'Failed to calculate custom weighting with uncertainty')
@@ -705,7 +706,7 @@ class Test(unittest.TestCase):
lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
grid_def = geometry.GridDefinition(lons, lats)
lons = numpy.asarray(lons, dtype='f4')
- lats = numpy.asarray(lats, dtype='f4')
+ lats = numpy.asarray(lats, dtype='f4')
swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
valid_input_index, valid_output_index, index_array, distance_array = \
kd_tree.get_neighbour_info(swath_def,
@@ -824,3 +825,6 @@ def suite():
mysuite.addTest(loader.loadTestsFromTestCase(Test))
return mysuite
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/pyresample/test/test_utils.py b/pyresample/test/test_utils.py
index 7f0409e..b43655c 100644
--- a/pyresample/test/test_utils.py
+++ b/pyresample/test/test_utils.py
@@ -221,6 +221,14 @@ class TestMisc(unittest.TestCase):
np.testing.assert_almost_equal(a, 6378137.)
np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)
+ 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_str_to_dict(proj_str)
+ proj_str2 = utils.proj4_dict_to_str(proj_dict)
+ proj_dict2 = utils.proj4_str_to_dict(proj_str2)
+ self.assertDictEqual(proj_dict, proj_dict2)
+
def suite():
"""The test suite.
diff --git a/pyresample/utils.py b/pyresample/utils.py
index 3351b54..7425481 100644
--- a/pyresample/utils.py
+++ b/pyresample/utils.py
@@ -105,11 +105,13 @@ def _read_yaml_area_file_content(area_file_name):
area_dict = {}
for area_file_obj in area_file_name:
if (isinstance(area_file_obj, (str, six.text_type)) and
- os.path.isfile(area_file_obj)):
- # filename
- area_file_obj = open(area_file_obj)
- tmp_dict = yaml.load(area_file_obj)
+ os.path.isfile(area_file_obj)):
+ with open(area_file_obj) as area_file_obj:
+ tmp_dict = yaml.load(area_file_obj)
+ else:
+ tmp_dict = yaml.load(area_file_obj)
area_dict = recursive_dict_update(area_dict, tmp_dict)
+
return area_dict
@@ -174,10 +176,10 @@ def _read_legacy_area_file_lines(area_file_name):
continue
elif isinstance(area_file_obj, (str, six.text_type)):
# filename
- area_file_obj = open(area_file_obj, 'r')
+ with open(area_file_obj, 'r') as area_file_obj:
- for line in area_file_obj.readlines():
- yield line
+ for line in area_file_obj.readlines():
+ yield line
def _parse_legacy_area_file(area_file_name, *regions):
@@ -281,11 +283,12 @@ def get_area_def(area_id, area_name, proj_id, proj4_args, x_size, y_size,
"""
proj_dict = _get_proj4_args(proj4_args)
- return pr.geometry.AreaDefinition(area_id, area_name, proj_id, proj_dict, x_size,
- y_size, area_extent)
+ return pr.geometry.AreaDefinition(area_id, area_name, proj_id, proj_dict,
+ x_size, y_size, area_extent)
-def generate_quick_linesample_arrays(source_area_def, target_area_def, nprocs=1):
+def generate_quick_linesample_arrays(source_area_def, target_area_def,
+ nprocs=1):
"""Generate linesample arrays for quick grid resampling
Parameters
@@ -315,8 +318,10 @@ def generate_quick_linesample_arrays(source_area_def, target_area_def, nprocs=1)
return source_pixel_y, source_pixel_x
-def generate_nearest_neighbour_linesample_arrays(source_area_def, target_area_def,
- radius_of_influence, nprocs=1):
+def generate_nearest_neighbour_linesample_arrays(source_area_def,
+ target_area_def,
+ radius_of_influence,
+ nprocs=1):
"""Generate linesample arrays for nearest neighbour grid resampling
Parameters
@@ -406,10 +411,31 @@ def proj4_str_to_dict(proj4_str):
Note: Key only parameters will be assigned a value of `True`.
"""
- pairs = (x.split('=', 1) for x in proj4_str.split(" "))
+ pairs = (x.split('=', 1) for x in proj4_str.replace('+', '').split(" "))
return dict((x[0], (x[1] if len(x) == 2 else True)) for x in pairs)
+def proj4_dict_to_str(proj4_dict, sort=False):
+ """Convert a dictionary of PROJ.4 parameters to a valid PROJ.4 string"""
+ keys = proj4_dict.keys()
+ if sort:
+ keys = sorted(keys)
+ params = []
+ for key in keys:
+ val = proj4_dict[key]
+ key = str(key) if key.startswith('+') else '+' + str(key)
+ if str(val) in ['True', 'False']:
+ # could be string or boolean object
+ val = ''
+ if val:
+ param = '{}={}'.format(key, val)
+ else:
+ # example "+no_defs"
+ param = key
+ params.append(param)
+ return ' '.join(params)
+
+
def proj4_radius_parameters(proj4_dict):
"""Calculate 'a' and 'b' radius parameters.
@@ -423,22 +449,30 @@ def proj4_radius_parameters(proj4_dict):
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']
+
+ from pyproj import Geod
+
+ if 'ellps' in new_info:
+ geod = Geod(**new_info)
+ new_info['a'] = geod.a
+ new_info['b'] = geod.b
+ elif 'a' not in new_info or 'b' not in new_info:
+
+ if 'rf' in new_info and 'f' not in new_info:
+ new_info['f'] = 1. / float(new_info['rf'])
+
+ if 'a' in new_info and 'f' in new_info:
+ new_info['b'] = float(new_info['a']) * (1 - float(new_info['f']))
+ elif 'b' in new_info and 'f' in new_info:
+ new_info['a'] = float(new_info['b']) / (1 - float(new_info['f']))
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'])
+ geod = Geod(**{'ellps': 'WGS84'})
+ new_info['a'] = geod.a
+ new_info['b'] = geod.b
+
+ return float(new_info['a']), float(new_info['b'])
def _downcast_index_array(index_array, size):
diff --git a/pyresample/version.py b/pyresample/version.py
index 6f605a3..3e0309b 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.1'
+__version__ = '1.7.0'
--
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