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 <[email protected]> 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 _______________________________________________ Pkg-grass-devel mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-grass-devel
