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 52e36972e7a752a5aaae937fde1f82add1014644 Author: Antonio Valentino <[email protected]> Date: Fri Dec 22 12:34:15 2017 +0100 New upstream version 1.7.1 --- .bumpversion.cfg | 2 +- .github/ISSUE_TEMPLATE.md | 18 ++++ .github/PULL_REQUEST_TEMPLATE.md | 7 ++ .travis.yml | 3 + changelog.rst | 97 ++++++++++++++++++++- pyresample/geometry.py | 183 +++++++++++++++++++++++++++++---------- pyresample/grid.py | 6 +- pyresample/image.py | 2 + pyresample/kd_tree.py | 27 +++--- pyresample/test/test_bilinear.py | 99 +++++++++++---------- pyresample/test/test_geometry.py | 128 +++++++++++++++++++++++++++ pyresample/test/test_grid.py | 16 ++++ pyresample/test/test_kd_tree.py | 63 +++++++++----- pyresample/test/test_rotation.py | 97 +++++++++++++++++++++ pyresample/utils.py | 21 +++-- pyresample/version.py | 2 +- 16 files changed, 636 insertions(+), 135 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 064551a..e1d2934 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.7.0 +current_version = 1.7.1 commit = True tag = True diff --git a/.github/ISSUE_TEMPLATE.md b/.github/ISSUE_TEMPLATE.md new file mode 100644 index 0000000..18c258a --- /dev/null +++ b/.github/ISSUE_TEMPLATE.md @@ -0,0 +1,18 @@ +#### Code Sample, a minimal, complete, and verifiable piece of code + +```python +# Your code here + +``` +#### Problem description + +[this should also explain **why** the current behaviour is a problem and why the +expected output is a better solution.] + +#### Expected Output + +#### Actual Result, Traceback if applicable + +#### Versions of Python, package at hand and relevant dependencies + +Thank you for reporting an issue ! diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md new file mode 100644 index 0000000..670a685 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -0,0 +1,7 @@ +Please make the PR against the `develop` branch. + + - [ ] Closes #xxxx (remove if there is no corresponding issue, which should only be the case for minor changes) + - [ ] Tests added (for all bug fixes or enhancements) + - [ ] Tests passed (for all non-documentation changes) + - [ ] Passes ``git diff origin/develop **/*py | flake8 --diff`` (remove if you did not edit any Python files) + - [ ] Fully documented (remove if this change should not be visible to users, e.g., if it is an internal clean-up, or if this is part of a larger project that will be documented later) diff --git a/.travis.yml b/.travis.yml index f7c6604..f3e5fae 100644 --- a/.travis.yml +++ b/.travis.yml @@ -21,6 +21,9 @@ install: - 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 toolz +- pip install dask +- pip install xarray - pip install -e ".[pykdtree,quicklook]" - pip install coveralls script: diff --git a/changelog.rst b/changelog.rst index 090290a..2560491 100644 --- a/changelog.rst +++ b/changelog.rst @@ -2,6 +2,92 @@ Changelog ========= +v1.7.1 (2017-12-21) +------------------- +- update changelog. [davidh-ssec] +- Bump version: 1.7.0 → 1.7.1. [davidh-ssec] +- Merge pull request #88 from pytroll/bugfix-masked-target. [David + Hoese] + + Fix kdtree when target lons/lats are masked arrays +- Add test for masked valid_output_index fix. [davidh-ssec] +- Move bilinear test setup to a special method. [davidh-ssec] +- Fix kdtree when target lons/lats are masked arrays. [davidh-ssec] +- Merge pull request #89 from Funkensieper/fix-masks-in-get-resampled- + image. [David Hoese] + + Fix masks in grid.get_resampled_image +- Add test for mask preservation. [Stephan Finkensieper] +- Distinguish between ndarrays and masked arrays. [Stephan Finkensieper] +- Fix masks in grid.get_resampled_image. [Stephan Finkensieper] + + Use numpy.ma version of row_stack to prevent loosing the mask of + large images (rows > cut_off) + +- Add github templates. [Martin Raspaud] +- Merge pull request #84 from pytroll/feature-add-hash. [Martin Raspaud] + + Add hash method to AreaDefinition and SwathDefinition +- Fix dask array not being hashable in py3.x. [Martin Raspaud] +- Use identity checking instead of equality. [Martin Raspaud] +- Do not has the mask if it's empty. [Martin Raspaud] +- Bugfix geometry test. [Martin Raspaud] +- Replace hash value checks with type checks. [Martin Raspaud] + + The value can be different depending on the python version apparently. +- Add dask and xarray for testing on travis. [Martin Raspaud] +- Fix case of missing xarray dependency in the tests. [Martin Raspaud] +- Add __hash__ for SwathDefinitions, along with some unittests. [Martin + Raspaud] +- Add hash method to AreaDefinition. [davidh-ssec] + + Removes annoying log message when xarray/dask is missing + +- Merge branch 'feature-xarray-improvements' into develop. [Martin + Raspaud] + + Conflicts: + pyresample/geometry.py + +- Type coords to np.float. [Martin Raspaud] +- Add support for fill_value in nn search. [Martin Raspaud] +- Change the get_lonlats_dask interface to return a tuple. [Martin + Raspaud] +- Fix masking bad latitude values. [davidh-ssec] +- Fix consistency with numpy arrays. [davidh-ssec] +- Allow xarrays internally in geometry objects. [davidh-ssec] +- Merge remote-tracking branch 'origin/develop' into develop. [davidh- + ssec] + + # Conflicts: + # .travis.yml + +- Fix proj4 dict to string against recent changes to str to dict funcs. + [davidh-ssec] +- Change appveyor python 3.5 environments to python 3.6. [davidh-ssec] + + Also removes slack notification webhook which is no longer the + recommended way to post to slack from appveyor. + +- Exclude buggy version of matplotlib in travis tests. [davidh-ssec] +- Fix proj4 dict conversion test. [davidh-ssec] +- Use more descriptive variable names. [davidh-ssec] +- Add proj4_dict_to_str utility function. [davidh-ssec] + + Includes fixes for dynamic area definitions proj_id and + small performance improvement for projection coordinate generation + +- Merge pull request #83 from loreclem/master. [Martin Raspaud] + + Added ROTATION in an area definition +- Bugfix in get_area_def. [lorenzo clementi] +- Unit test for rotation. [lorenzo clementi] +- Removed unused parameter. [lorenzo clementi] +- Now working also with yaml. [lorenzo clementi] +- Code improvements. [lorenzo clementi] +- Added ROTATION in an area definition. [lorenzo clementi] + + v1.7.0 (2017-10-13) ------------------- - update changelog. [Martin Raspaud] @@ -526,6 +612,7 @@ v1.2.2 (2016-06-21) Without this, the compilation of the ewa extension crashes. + v1.2.1 (2016-06-21) ------------------- - update changelog. [Martin Raspaud] @@ -681,9 +768,11 @@ v1.2.0 (2016-06-17) - Make kd_tree test work on older numpy version. [Martin Raspaud] VisibleDeprecationWarning is not available in numpy <1.9. + - Adapt to newest pykdtree version. [Martin Raspaud] The kdtree object's attribute `data_pts` has been renamed to `data`. + - Run tests on python 3.5 in travis also. [Martin Raspaud] @@ -695,6 +784,7 @@ v1.1.6 (2016-02-25) A previous commit was looking for a 'data_pts' attribute in the kdtree object, which is available in pykdtree, but not scipy. + - Merge pull request #32 from mitkin/master. [Martin Raspaud] [tests] Skip deprecation warnings in test_gauss_multi_uncert @@ -706,6 +796,7 @@ v1.1.6 (2016-02-25) The latest matplotlib (1.5) doesn't support python 2.6 and 3.3. This patch chooses the right matplotlib version to install depending on the python version at hand. + - Skip deprecation warnings. [Mikhail Itkin] Catch the rest of the warnings. Check if there is only one, and @@ -747,6 +838,7 @@ Other - Bugfix to address a numpy DeprecationWarning. [Martin Raspaud] Numpy won't take non-integer indices soon, so make index an int. + - Merge branch 'release-1.1.3' [Martin Raspaud] - Merge branch 'licence-lgpl' into pre-master. [Martin Raspaud] - Switch to lgplv3, and bump up version number. [Martin Raspaud] @@ -968,7 +1060,7 @@ Other - Set svn:mime-type. [StorPipfugl] - Corrected doc errors. [StorPipfugl] - Removed dist dir. [StorPipfugl] -- No commit message. [StorPipfugl] +- [StorPipfugl] - Updated documentation. New release. [StorPipfugl] - Started updating docstrings. [StorPipfugl] - Restructured API. [StorPipfugl] @@ -981,9 +1073,8 @@ Other - Removed unneeded function. [StorPipfugl] - Mime types set. [StorPipfugl] - Mime types set. [StorPipfugl] -- No commit message. [StorPipfugl] +- [StorPipfugl] - Moved to Google Code under GPLv3 license. [StorPipfugl] - moved to Google Code. [StorPipfugl] - diff --git a/pyresample/geometry.py b/pyresample/geometry.py index 988c310..9044f0f 100644 --- a/pyresample/geometry.py +++ b/pyresample/geometry.py @@ -25,6 +25,7 @@ import warnings from collections import OrderedDict from logging import getLogger +import hashlib import numpy as np import yaml @@ -32,6 +33,11 @@ from pyproj import Geod, Proj from pyresample import _spatial_mp, utils +try: + from xarray import DataArray +except ImportError: + DataArray = np.ndarray + logger = getLogger(__name__) @@ -63,31 +69,37 @@ class BaseDefinition(object): if type(lons) != type(lats): raise TypeError('lons and lats must be of same type') elif lons is not None: - lons = np.asanyarray(lons) - lats = np.asanyarray(lats) + if not isinstance(lons, (np.ndarray, DataArray)): + lons = np.asanyarray(lons) + lats = np.asanyarray(lats) if lons.shape != lats.shape: raise ValueError('lons and lats must have same shape') self.nprocs = nprocs # check the latitutes - if lats is not None and ((lats.min() < -90. or lats.max() > +90.)): - # throw exception - raise ValueError( - 'Some latitudes are outside the [-90.;+90] validity range') - else: - self.lats = lats + if lats is not None: + if isinstance(lats, np.ndarray) and (lats.min() < -90. or lats.max() > 90.): + raise ValueError( + 'Some latitudes are outside the [-90.;+90] validity range') + elif not isinstance(lats, np.ndarray): + # assume we have to mask an xarray + lats = lats.where((lats >= -90.) & (lats <= 90.)) + self.lats = lats # check the longitudes - if lons is not None and ((lons.min() < -180. or lons.max() >= +180.)): - # issue warning - warnings.warn('All geometry objects expect longitudes in the [-180:+180[ range. ' + - 'We will now automatically wrap your longitudes into [-180:+180[, and continue. ' + - 'To avoid this warning next time, use routine utils.wrap_longitudes().') - # wrap longitudes to [-180;+180[ - self.lons = utils.wrap_longitudes(lons) - else: - self.lons = lons + if lons is not None: + if isinstance(lons, np.ndarray) and (lons.min() < -180. or lons.max() >= +180.): + # issue warning + warnings.warn('All geometry objects expect longitudes in the [-180:+180[ range. ' + + 'We will now automatically wrap your longitudes into [-180:+180[, and continue. ' + + 'To avoid this warning next time, use routine utils.wrap_longitudes().') + # assume we have to mask an xarray + # wrap longitudes to [-180;+180[ + lons = utils.wrap_longitudes(lons) + elif not isinstance(lons, np.ndarray): + lons = utils.wrap_longitudes(lons) + self.lons = lons self.ndim = None self.cartesian_coords = None @@ -357,8 +369,9 @@ class CoordinateDefinition(BaseDefinition): """Base class for geometry definitions defined by lons and lats only""" def __init__(self, lons, lats, nprocs=1): - lons = np.asanyarray(lons) - lats = np.asanyarray(lats) + if not isinstance(lons, (np.ndarray, DataArray)): + lons = np.asanyarray(lons) + lats = np.asanyarray(lats) super(CoordinateDefinition, self).__init__(lons, lats, nprocs) if lons.shape == lats.shape and lons.dtype == lats.dtype: self.shape = lons.shape @@ -429,9 +442,23 @@ class GridDefinition(CoordinateDefinition): raise ValueError('2 dimensional lon lat grid expected') -class SwathDefinition(CoordinateDefinition): +def get_array_hashable(arr): + """Compute a hashable form of the array `arr`. + + Works with numpy arrays, dask.array.Array, and xarray.DataArray. + """ + # look for precomputed value + if isinstance(arr, DataArray) and np.ndarray is not DataArray: + return arr.attrs.get('hash', get_array_hashable(arr.data)) + else: + try: + return arr.name.encode('utf-8') # dask array + except AttributeError: + return arr.view(np.uint8) # np array - """Swath defined by lons and lats + +class SwathDefinition(CoordinateDefinition): + """Swath defined by lons and lats. Parameters ---------- @@ -454,24 +481,49 @@ class SwathDefinition(CoordinateDefinition): Swath lats cartesian_coords : object Swath cartesian coordinates + """ def __init__(self, lons, lats, nprocs=1): - lons = np.asanyarray(lons) - lats = np.asanyarray(lats) + if not isinstance(lons, (np.ndarray, DataArray)): + lons = np.asanyarray(lons) + lats = np.asanyarray(lats) super(SwathDefinition, self).__init__(lons, lats, nprocs) if lons.shape != lats.shape: raise ValueError('lon and lat arrays must have same shape') elif lons.ndim > 2: raise ValueError('Only 1 and 2 dimensional swaths are allowed') + self.hash = None + + def __hash__(self): + """Compute the hash of this object.""" + if self.hash is None: + hasher = hashlib.sha1() + hasher.update(get_array_hashable(self.lons)) + hasher.update(get_array_hashable(self.lats)) + try: + if self.lons.mask is not np.bool_(False): + hasher.update(get_array_hashable(self.lons.mask)) + except AttributeError: + pass + self.hash = int(hasher.hexdigest(), 16) + + return self.hash + 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) + + if isinstance(lons.data, da.Array): + return lons.data, lats.data + else: + lons = da.from_array(np.asanyarray(lons), + chunks=blocksize * lons.ndim) + lats = da.from_array(np.asanyarray(lats), + chunks=blocksize * lats.ndim) + return lons, lats def _compute_omerc_parameters(self, ellipsoid): """Compute the oblique mercator projection bouding box parameters.""" @@ -538,13 +590,14 @@ class SwathDefinition(CoordinateDefinition): class DynamicAreaDefinition(object): """An AreaDefintion containing just a subset of the needed parameters. - The purpose of this class is to be able to adapt the area extent and size of - the area to a given set of longitudes and latitudes, such that e.g. polar - satellite granules can be resampled optimaly to a give projection.""" + The purpose of this class is to be able to adapt the area extent and size + of the area to a given set of longitudes and latitudes, such that e.g. + polar satellite granules can be resampled optimaly to a give projection. + """ def __init__(self, area_id=None, description=None, proj_dict=None, x_size=None, y_size=None, area_extent=None, - optimize_projection=False): + optimize_projection=False, rotation=None): """Initialize the DynamicAreaDefinition. area_id: @@ -559,6 +612,8 @@ class DynamicAreaDefinition(object): The area extent of the area. optimize_projection: Whether the projection parameters have to be optimized. + rotation: + Rotation in degrees (negative is cw) """ self.area_id = area_id self.description = description @@ -567,6 +622,7 @@ class DynamicAreaDefinition(object): self.y_size = y_size self.area_extent = area_extent self.optimize_projection = optimize_projection + self.rotation = rotation def compute_domain(self, corners, resolution=None, size=None): """Compute size and area_extent from corners and [size or resolution] @@ -597,7 +653,7 @@ class DynamicAreaDefinition(object): def freeze(self, lonslats=None, resolution=None, size=None, - proj_info=None): + proj_info=None, rotation=None): """Create an AreaDefintion from this area with help of some extra info. lonlats: @@ -608,6 +664,8 @@ class DynamicAreaDefinition(object): the size of the resulting area. proj_info: complementing parameters to the projection info. + rotation: + rotation in degrees (negative is cw) Resolution and Size parameters are ignored if the instance is created with the `optimize_projection` flag set to True. @@ -632,7 +690,7 @@ class DynamicAreaDefinition(object): return AreaDefinition(self.area_id, self.description, '', self.proj_dict, self.x_size, self.y_size, - self.area_extent) + self.area_extent, self.rotation) class AreaDefinition(BaseDefinition): @@ -653,6 +711,8 @@ class AreaDefinition(BaseDefinition): x dimension in number of pixels y_size : int y dimension in number of pixels + rotation: float + rotation in degrees (negative is cw) area_extent : list Area extent as a list (LL_x, LL_y, UR_x, UR_y) nprocs : int, optional @@ -676,6 +736,8 @@ class AreaDefinition(BaseDefinition): x dimension in number of pixels y_size : int y dimension in number of pixels + rotation: float + rotation in degrees (negative is cw) shape : tuple Corresponding array shape as (rows, cols) size : int @@ -711,7 +773,8 @@ class AreaDefinition(BaseDefinition): """ def __init__(self, area_id, name, proj_id, proj_dict, x_size, y_size, - area_extent, nprocs=1, lons=None, lats=None, dtype=np.float64): + area_extent, rotation=None, nprocs=1, lons=None, lats=None, + dtype=np.float64): if not isinstance(proj_dict, dict): raise TypeError('Wrong type for proj_dict: %s. Expected dict.' % type(proj_dict)) @@ -723,6 +786,10 @@ class AreaDefinition(BaseDefinition): self.x_size = int(x_size) self.y_size = int(y_size) self.shape = (y_size, x_size) + try: + self.rotation = float(rotation) + except TypeError: + self.rotation = 0 if lons is not None: if lons.shape != self.shape: raise ValueError('Shape of lon lat grid must match ' @@ -749,8 +816,8 @@ class AreaDefinition(BaseDefinition): float(area_extent[3]) - float(self.pixel_size_y) / 2) - # Pixel_offset defines the distance to projection center from origen (UL) - # of image in units of pixels. + # Pixel_offset defines the distance to projection center from origen + # (UL) of image in units of pixels. self.pixel_offset_x = -self.area_extent[0] / self.pixel_size_x self.pixel_offset_y = self.area_extent[3] / self.pixel_size_y @@ -810,6 +877,7 @@ class AreaDefinition(BaseDefinition): fmt += "\tPCS_DEF:\t{proj_str}\n" fmt += "\tXSIZE:\t{x_size}\n" fmt += "\tYSIZE:\t{y_size}\n" + fmt += "\tROTATION:\t{rotation}\n" fmt += "\tAREA_EXTENT: {area_extent}\n}};\n" area_def_str = fmt.format(name=self.name, area_id=self.area_id, proj_str=proj_str, x_size=self.x_size, @@ -823,7 +891,7 @@ class AreaDefinition(BaseDefinition): """Test for equality""" try: - return ((self.proj_dict == other.proj_dict) and + return ((self.proj_str == other.proj_str) and (self.shape == other.shape) and (np.allclose(self.area_extent, other.area_extent))) except AttributeError: @@ -834,6 +902,13 @@ class AreaDefinition(BaseDefinition): return not self.__eq__(other) + def __hash__(self): + return hash(( + self.proj_str, + self.shape, + self.area_extent + )) + def colrow2lonlat(self, cols, rows): """ Return longitudes and latitudes for the given image columns @@ -1007,6 +1082,13 @@ class AreaDefinition(BaseDefinition): Grids of area x- and y-coordinates in projection units """ + def do_rotation(xspan, yspan, rot_deg=0): + rot_rad = np.radians(rot_deg) + rot_mat = np.array([[np.cos(rot_rad), np.sin(rot_rad)], + [-np.sin(rot_rad), np.cos(rot_rad)]]) + x, y = np.meshgrid(xspan, yspan) + return np.einsum('ji, mni -> jmn', rot_mat, np.dstack([x, y])) + def get_val(val, sub_val, max): # Get value with substitution and wrapping if val is None: @@ -1081,9 +1163,15 @@ class AreaDefinition(BaseDefinition): cols = 1 # Calculate coordinates - 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) + 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] + if self.rotation != 0: + res = do_rotation(target_x, target_y, self.rotation) + target_x, target_y = res[0, :, :], res[1, :, :] + else: + target_x, target_y = np.meshgrid(target_x, target_y) if is_single_value: # Return single values @@ -1169,7 +1257,7 @@ class AreaDefinition(BaseDefinition): res = Array(dsk, name, shape=list(self.shape) + [2], chunks=(blocksize, blocksize, 2), dtype=dtype) - return res + return res[:, :, 0], res[:, :, 1] def get_lonlats(self, nprocs=None, data_slice=None, cache=False, dtype=None): """Returns lon and lat arrays of area. @@ -1363,16 +1451,19 @@ class StackedAreaDefinition(BaseDefinition): def get_lonlats_dask(self, blocksize=1000, dtype=None): """"Return lon and lat dask arrays of the area.""" import dask.array as da - llonslats = [] + llons = [] + llats = [] for definition in self.defs: - lonslats = definition.get_lonlats_dask(blocksize=blocksize, - dtype=dtype) + lons, lats = definition.get_lonlats_dask(blocksize=blocksize, + dtype=dtype) - llonslats.append(lonslats) + llons.append(lons) + llats.append(lats) - self.lonlats = da.concatenate(llonslats, axis=0) + self.lons = da.concatenate(llons, axis=0) + self.lats = da.concatenate(llats, axis=0) - return self.lonlats + return self.lons, self.lats def squeeze(self): """Generate a single AreaDefinition if possible.""" diff --git a/pyresample/grid.py b/pyresample/grid.py index a484e59..05afb1e 100644 --- a/pyresample/grid.py +++ b/pyresample/grid.py @@ -236,7 +236,11 @@ def get_resampled_image(target_area_def, source_area_def, source_image_data, # First iteration result = next_result else: - result = np.row_stack((result, next_result)) + if isinstance(next_result, np.ma.core.MaskedArray): + stack = np.ma.row_stack + else: + stack = np.row_stack + result = stack((result, next_result)) return result else: diff --git a/pyresample/image.py b/pyresample/image.py index 5254ac1..629b2d7 100644 --- a/pyresample/image.py +++ b/pyresample/image.py @@ -55,6 +55,8 @@ class ImageContainer(object): """ def __init__(self, image_data, geo_def, fill_value=0, nprocs=1): + if type(geo_def).__name__ == "DynamicAreaDefinition": + geo_def = geo_def.freeze() if not isinstance(image_data, (np.ndarray, np.ma.core.MaskedArray)): raise TypeError('image_data must be either an ndarray' ' or a masked array') diff --git a/pyresample/kd_tree.py b/pyresample/kd_tree.py index a838461..0f9f54c 100644 --- a/pyresample/kd_tree.py +++ b/pyresample/kd_tree.py @@ -36,7 +36,8 @@ try: from xarray import DataArray import dask.array as da except ImportError: - logger.info("XArray or dask unavailable, some functionality missing.") + DataArray = None + da = None if sys.version < '3': range = xrange @@ -466,6 +467,8 @@ def _get_valid_output_index(source_geo_def, target_geo_def, target_lons, # Combine reduced and legal values valid_output_index = (valid_output_index & valid_out) + if isinstance(valid_output_index, np.ma.MaskedArray): + valid_output_index = valid_output_index.filled(False) return valid_output_index @@ -914,6 +917,8 @@ class XArrayResamplerNN(object): Number of segments to use when resampling. If set to None an estimate will be calculated """ + if DataArray is None: + raise ImportError("Missing 'xarray' and 'dask' dependencies") self.valid_input_index = None self.valid_output_index = None @@ -960,7 +965,7 @@ class XArrayResamplerNN(object): raise EmptyResult('No valid data points in input data') # Build kd-tree on input - + input_coords = input_coords.astype(np.float) if kd_tree_name == 'pykdtree': resample_kdtree = KDTree(input_coords.compute()) else: @@ -1036,9 +1041,7 @@ class XArrayResamplerNN(object): 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] + source_lons, source_lats = self.source_geo_def.get_lonlats_dask() valid_input_index = ((source_lons >= -180) & (source_lons <= 180) & (source_lats <= 90) & (source_lats >= -90)) @@ -1053,12 +1056,14 @@ class XArrayResamplerNN(object): valid_output_index, index_array, distance_array = \ _create_empty_info(self.source_geo_def, self.target_geo_def, self.neighbours) + 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) - target_lonlats = self.target_geo_def.get_lonlats_dask() - target_lons = target_lonlats[:, :, 0] - target_lats = target_lonlats[:, :, 1] + target_lons, target_lats = self.target_geo_def.get_lonlats_dask() valid_output_index = ((target_lons >= -180) & (target_lons <= 180) & (target_lats <= 90) & (target_lats >= -90)) @@ -1074,7 +1079,7 @@ class XArrayResamplerNN(object): return valid_input_index, valid_output_index, index_array, distance_array - def get_sample_from_neighbour_info(self, data): + def get_sample_from_neighbour_info(self, data, fill_value=np.nan): # flatten x and y in the source array @@ -1124,9 +1129,9 @@ class XArrayResamplerNN(object): 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 + result[index_mask.ravel()] = fill_value #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 = np.full(target_shape[0], fill_value) target_data_line[valid_targets] = result target_lines.append(target_data_line[:, np.newaxis]) diff --git a/pyresample/test/test_bilinear.py b/pyresample/test/test_bilinear.py index 7a0eb1c..f5a9051 100644 --- a/pyresample/test/test_bilinear.py +++ b/pyresample/test/test_bilinear.py @@ -9,52 +9,59 @@ from pyresample import geometry, utils, kd_tree class Test(unittest.TestCase): - pts_irregular = (np.array([[-1., 1.], ]), - np.array([[1., 2.], ]), - np.array([[-2., -1.], ]), - np.array([[2., -4.], ])) - pts_vert_parallel = (np.array([[-1., 1.], ]), - np.array([[1., 2.], ]), - np.array([[-1., -1.], ]), - np.array([[1., -2.], ])) - pts_both_parallel = (np.array([[-1., 1.], ]), - np.array([[1., 1.], ]), - np.array([[-1., -1.], ]), - np.array([[1., -1.], ])) - - # Area definition with four pixels - target_def = geometry.AreaDefinition('areaD', - 'Europe (3km, HRV, VTC)', - 'areaD', - {'a': '6378144.0', - 'b': '6356759.0', - 'lat_0': '50.00', - 'lat_ts': '50.00', - 'lon_0': '8.00', - 'proj': 'stere'}, - 4, 4, - [-1370912.72, - -909968.64000000001, - 1029087.28, - 1490031.3600000001]) - - # Input data around the target pixel at 0.63388324, 55.08234642, - in_shape = (100, 100) - data1 = np.ones((in_shape[0], in_shape[1])) - data2 = 2. * data1 - lons, lats = np.meshgrid(np.linspace(-5., 5., num=in_shape[0]), - np.linspace(50., 60., num=in_shape[1])) - swath_def = geometry.SwathDefinition(lons=lons, lats=lats) - - radius = 50e3 - neighbours = 32 - input_idxs, output_idxs, idx_ref, dists = \ - kd_tree.get_neighbour_info(swath_def, target_def, - radius, neighbours=neighbours, - nprocs=1) - input_size = input_idxs.sum() - index_mask = (idx_ref == input_size) - idx_ref = np.where(index_mask, 0, idx_ref) + @classmethod + def setUpClass(cls): + cls.pts_irregular = (np.array([[-1., 1.], ]), + np.array([[1., 2.], ]), + np.array([[-2., -1.], ]), + np.array([[2., -4.], ])) + cls.pts_vert_parallel = (np.array([[-1., 1.], ]), + np.array([[1., 2.], ]), + np.array([[-1., -1.], ]), + np.array([[1., -2.], ])) + cls.pts_both_parallel = (np.array([[-1., 1.], ]), + np.array([[1., 1.], ]), + np.array([[-1., -1.], ]), + np.array([[1., -1.], ])) + + # Area definition with four pixels + target_def = geometry.AreaDefinition('areaD', + 'Europe (3km, HRV, VTC)', + 'areaD', + {'a': '6378144.0', + 'b': '6356759.0', + 'lat_0': '50.00', + 'lat_ts': '50.00', + 'lon_0': '8.00', + 'proj': 'stere'}, + 4, 4, + [-1370912.72, + -909968.64000000001, + 1029087.28, + 1490031.3600000001]) + + # Input data around the target pixel at 0.63388324, 55.08234642, + in_shape = (100, 100) + cls.data1 = np.ones((in_shape[0], in_shape[1])) + cls.data2 = 2. * cls.data1 + lons, lats = np.meshgrid(np.linspace(-5., 5., num=in_shape[0]), + np.linspace(50., 60., num=in_shape[1])) + cls.swath_def = geometry.SwathDefinition(lons=lons, lats=lats) + + radius = 50e3 + cls.neighbours = 32 + input_idxs, output_idxs, idx_ref, dists = \ + kd_tree.get_neighbour_info(cls.swath_def, target_def, + radius, neighbours=cls.neighbours, + nprocs=1) + input_size = input_idxs.sum() + index_mask = (idx_ref == input_size) + idx_ref = np.where(index_mask, 0, idx_ref) + + cls.input_idxs = input_idxs + cls.target_def = target_def + cls.idx_ref = idx_ref + def test_calc_abc(self): # No np.nan inputs diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py index 73a3050..4aec478 100644 --- a/pyresample/test/test_geometry.py +++ b/pyresample/test/test_geometry.py @@ -165,6 +165,134 @@ class Test(unittest.TestCase): "BaseDefinition did not maintain dtype of longitudes (in:%s out:%s)" % (lons2_ints.dtype, lons.dtype,)) + def test_area_hash(self): + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + {'a': '6378144.0', + 'b': '6356759.0', + 'lat_0': '50.00', + 'lat_ts': '50.00', + 'lon_0': '8.00', + 'proj': 'stere'}, + 800, + 800, + [-1370912.72, + -909968.64000000001, + 1029087.28, + 1490031.3600000001]) + + self.assertIsInstance(hash(area_def), int) + + area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', + {'a': '6378144.0', + 'b': '6356759.0', + 'lat_ts': '50.00', + 'lon_0': '8.00', + 'lat_0': '50.00', + 'proj': 'stere'}, + 800, + 800, + [-1370912.72, + -909968.64000000001, + 1029087.28, + 1490031.3600000001]) + + self.assertIsInstance(hash(area_def), int) + + area_def = geometry.AreaDefinition('New area', 'Europe', 'areaD', + {'a': '6378144.0', + 'b': '6356759.0', + 'lat_ts': '50.00', + 'lon_0': '8.00', + 'lat_0': '50.00', + 'proj': 'stere'}, + 800, + 800, + [-1370912.72, + -909968.64000000001, + 1029087.28, + 1490031.3600000001]) + + self.assertIsInstance(hash(area_def), int) + + def test_get_array_hashable(self): + arr = np.array([1.2, 1.3, 1.4, 1.5]) + + np.testing.assert_allclose(np.array([ 51, 51, 51, 51, 51, 51, 243, + 63, 205, 204, 204, 204, 204, + 204, 244, 63, 102, 102, 102, 102, + 102, 102, 246, 63, 0, 0, + 0, 0, 0, 0, 248, 63], + dtype=np.uint8), + geometry.get_array_hashable(arr)) + + try: + import xarray as xr + except ImportError: + pass + else: + xrarr = xr.DataArray(arr) + np.testing.assert_allclose(np.array([ 51, 51, 51, 51, 51, 51, 243, + 63, 205, 204, 204, 204, 204, + 204, 244, 63, 102, 102, 102, 102, + 102, 102, 246, 63, 0, 0, + 0, 0, 0, 0, 248, 63], + dtype=np.uint8), + geometry.get_array_hashable(arr)) + + xrarr.attrs['hash'] = 42 + self.assertEqual(geometry.get_array_hashable(xrarr), + xrarr.attrs['hash']) + + + def test_swath_hash(self): + lons = np.array([1.2, 1.3, 1.4, 1.5]) + lats = np.array([65.9, 65.86, 65.82, 65.78]) + swath_def = geometry.SwathDefinition(lons, lats) + + self.assertIsInstance(hash(swath_def), int) + + try: + import dask.array as da + except ImportError: + print("Not testing with dask arrays") + else: + dalons = da.from_array(lons, chunks=1000) + dalats = da.from_array(lats, chunks=1000) + swath_def = geometry.SwathDefinition(dalons, dalats) + + self.assertIsInstance(hash(swath_def), int) + + try: + import xarray as xr + except ImportError: + print("Not testing with xarray") + else: + xrlons = xr.DataArray(lons) + xrlats = xr.DataArray(lats) + swath_def = geometry.SwathDefinition(xrlons, xrlats) + + self.assertIsInstance(hash(swath_def), int) + + try: + import xarray as xr + import dask.array as da + except ImportError: + print("Not testing with xarrays and dask arrays") + else: + xrlons = xr.DataArray(da.from_array(lons, chunks=1000)) + xrlats = xr.DataArray(da.from_array(lats, chunks=1000)) + swath_def = geometry.SwathDefinition(xrlons, xrlats) + + self.assertIsInstance(hash(swath_def), int) + + + lons = np.ma.array([1.2, 1.3, 1.4, 1.5]) + lats = np.ma.array([65.9, 65.86, 65.82, 65.78]) + swath_def = geometry.SwathDefinition(lons, lats) + + self.assertIsInstance(hash(swath_def), int) + + def test_area_equal(self): area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', {'a': '6378144.0', diff --git a/pyresample/test/test_grid.py b/pyresample/test/test_grid.py index 9a7f53d..9abd5f4 100644 --- a/pyresample/test/test_grid.py +++ b/pyresample/test/test_grid.py @@ -152,6 +152,22 @@ class Test(unittest.TestCase): self.assertAlmostEqual( cross_sum, expected, msg='Resampling of image failed') + def test_resampled_image_masked(self): + # Generate test image with masked elements + data = np.ma.ones(self.msg_area.shape) + data.mask = np.zeros(data.shape) + data.mask[253:400, 1970:2211] = 1 + + # Resample image using multiple segments + target_def = self.area_def + source_def = self.msg_area + res = grid.get_resampled_image( + target_def, source_def, data, segments=4, fill_value=None) + + # Make sure the mask has been preserved + self.assertGreater(res.mask.sum(), 0, + msg='Resampling did not preserve the mask') + @tmp def test_generate_linesample(self): data = np.fromfunction(lambda y, x: y * x * 10 ** -6, (3712, 3712)) diff --git a/pyresample/test/test_kd_tree.py b/pyresample/test/test_kd_tree.py index 1ec401a..d1cadd2 100644 --- a/pyresample/test/test_kd_tree.py +++ b/pyresample/test/test_kd_tree.py @@ -16,26 +16,30 @@ else: class Test(unittest.TestCase): - area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', - {'a': '6378144.0', - 'b': '6356759.0', - 'lat_0': '50.00', - 'lat_ts': '50.00', - 'lon_0': '8.00', - 'proj': 'stere'}, - 800, - 800, - [-1370912.72, - -909968.64000000001, - 1029087.28, - 1490031.3600000001]) - - tdata = numpy.array([1, 2, 3]) - tlons = numpy.array([11.280789, 12.649354, 12.080402]) - tlats = numpy.array([56.011037, 55.629675, 55.641535]) - tswath = geometry.SwathDefinition(lons=tlons, lats=tlats) - tgrid = geometry.CoordinateDefinition(lons=numpy.array([12.562036]), - lats=numpy.array([55.715613])) + @classmethod + def setUpClass(cls): + cls.area_def = geometry.AreaDefinition('areaD', + 'Europe (3km, HRV, VTC)', + 'areaD', + {'a': '6378144.0', + 'b': '6356759.0', + 'lat_0': '50.00', + 'lat_ts': '50.00', + 'lon_0': '8.00', + 'proj': 'stere'}, + 800, + 800, + [-1370912.72, + -909968.64000000001, + 1029087.28, + 1490031.3600000001]) + + cls.tdata = numpy.array([1, 2, 3]) + cls.tlons = numpy.array([11.280789, 12.649354, 12.080402]) + cls.tlats = numpy.array([56.011037, 55.629675, 55.641535]) + cls.tswath = geometry.SwathDefinition(lons=cls.tlons, lats=cls.tlats) + cls.tgrid = geometry.CoordinateDefinition( + lons=numpy.array([12.562036]), lats=numpy.array([55.715613])) def test_nearest_base(self): res = kd_tree.resample_nearest(self.tswath, @@ -123,6 +127,25 @@ class Test(unittest.TestCase): self.assertEqual(cross_sum, expected, msg='Swath resampling nearest failed') + def test_nearest_masked_swath_target(self): + """Test that a masked array works as a target.""" + 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)) + mask = numpy.ones_like(lons, dtype=numpy.bool) + mask[::2, ::2] = False + swath_def = geometry.SwathDefinition( + lons=numpy.ma.masked_array(lons, mask=mask), # numpy.ones_like(lons, dtype=numpy.bool)), + lats=numpy.ma.masked_array(lats, mask=False) + ) + res = kd_tree.resample_nearest(swath_def, data.ravel(), + swath_def, 50000, segments=3) + cross_sum = res.sum() + # expected = 12716 # if masks aren't respected + expected = 12000 + self.assertEqual(cross_sum, expected, + msg='Swath resampling masked nearest failed') + def test_nearest_1d(self): data = numpy.fromfunction(lambda x, y: x * y, (800, 800)) lons = numpy.fromfunction(lambda x: 3 + x / 100., (500,)) diff --git a/pyresample/test/test_rotation.py b/pyresample/test/test_rotation.py new file mode 100644 index 0000000..0a8c4db --- /dev/null +++ b/pyresample/test/test_rotation.py @@ -0,0 +1,97 @@ +import unittest, os +from pyresample.utils import load_area + +class TestProjRotation(unittest.TestCase): + + def test_rotation_legacy (self): + legacyDef = """REGION: regionB { + NAME: regionB + PCS_ID: regionB + PCS_DEF: proj=merc, lon_0=-34, k=1, x_0=0, y_0=0, a=6378137, b=6378137 + XSIZE: 800 + YSIZE: 548 + ROTATION: -45 + AREA_EXTENT: (-7761424.714818418, -4861746.639279127, 11136477.43264252, 8236799.845095873) + };""" + flegacy = "/tmp/TestProjRotation_test_rotation_legacy.txt" + f = open(flegacy,"w") + f.write(legacyDef) + f.close() + test_area = load_area(flegacy, 'regionB') + self.assertEqual(test_area.rotation, -45) + os.remove(flegacy) + + def test_rotation_yaml (self): + yamlDef = """regionB: + description: regionB + projection: + a: 6378137.0 + b: 6378137.0 + lon_0: -34 + proj: merc + x_0: 0 + y_0: 0 + k_0: 1 + shape: + height: 548 + width: 800 + rotation: -45 + area_extent: + lower_left_xy: [-7761424.714818418, -4861746.639279127] + upper_right_xy: [11136477.43264252, 8236799.845095873] + units: m""" + fyaml = "/tmp/TestProjRotation_test_rotation_yaml.txt" + f = open(fyaml,"w") + f.write(yamlDef) + f.close() + test_area = load_area(fyaml, 'regionB') + self.assertEqual(test_area.rotation, -45) + os.remove(fyaml) + + def test_norotation_legacy (self): + legacyDef = """REGION: regionB { + NAME: regionB + PCS_ID: regionB + PCS_DEF: proj=merc, lon_0=-34, k=1, x_0=0, y_0=0, a=6378137, b=6378137 + XSIZE: 800 + YSIZE: 548 + AREA_EXTENT: (-7761424.714818418, -4861746.639279127, 11136477.43264252, 8236799.845095873) + };""" + flegacy = "/tmp/TestProjRotation_test_rotation_legacy.txt" + f = open(flegacy,"w") + f.write(legacyDef) + f.close() + test_area = load_area(flegacy, 'regionB') + self.assertEqual(test_area.rotation, 0) + os.remove(flegacy) + + def test_norotation_yaml (self): + yamlDef = """regionB: + description: regionB + projection: + a: 6378137.0 + b: 6378137.0 + lon_0: -34 + proj: merc + x_0: 0 + y_0: 0 + k_0: 1 + shape: + height: 548 + width: 800 + area_extent: + lower_left_xy: [-7761424.714818418, -4861746.639279127] + upper_right_xy: [11136477.43264252, 8236799.845095873] + units: m""" + fyaml = "/tmp/TestProjRotation_test_rotation_yaml.txt" + f = open(fyaml,"w") + f.write(yamlDef) + f.close() + test_area = load_area(fyaml, 'regionB') + self.assertEqual(test_area.rotation, 0) + os.remove(fyaml) + + + +if __name__ == '__main__': + unittest.main() diff --git a/pyresample/utils.py b/pyresample/utils.py index 7425481..c97ae5a 100644 --- a/pyresample/utils.py +++ b/pyresample/utils.py @@ -150,10 +150,15 @@ def _parse_yaml_area_file(area_file_name, *regions): params['area_extent']['upper_right_xy']) except KeyError: area_extent = None + try: + rotation = params['rotation'] + except KeyError: + rotation = 0 area = pr.geometry.DynamicAreaDefinition(area_name, description, projection, xsize, ysize, area_extent, - optimize_projection) + optimize_projection, + rotation) try: area = area.freeze() except (TypeError, AttributeError): @@ -241,6 +246,10 @@ def _create_area(area_id, area_content): config['XSIZE'] = int(config['XSIZE']) config['YSIZE'] = int(config['YSIZE']) + if 'ROTATION' in config.keys(): + config['ROTATION'] = float(config['ROTATION']) + else: + config['ROTATION'] = 0 config['AREA_EXTENT'][0] = config['AREA_EXTENT'][0].replace('(', '') config['AREA_EXTENT'][3] = config['AREA_EXTENT'][3].replace(')', '') @@ -248,15 +257,14 @@ def _create_area(area_id, area_content): config['AREA_EXTENT'][i] = float(val) config['PCS_DEF'] = _get_proj4_args(config['PCS_DEF']) - return pr.geometry.AreaDefinition(config['REGION'], config['NAME'], config['PCS_ID'], config['PCS_DEF'], config['XSIZE'], config['YSIZE'], - config['AREA_EXTENT']) + config['AREA_EXTENT'], config['ROTATION']) def get_area_def(area_id, area_name, proj_id, proj4_args, x_size, y_size, - area_extent): + area_extent, rotation=0): """Construct AreaDefinition object from arguments Parameters @@ -273,6 +281,8 @@ def get_area_def(area_id, area_name, proj_id, proj4_args, x_size, y_size, Number of pixel in x dimension y_size : int Number of pixel in y dimension + rotation: float + Rotation in degrees (negative is cw) area_extent : list Area extent as a list of ints (LL_x, LL_y, UR_x, UR_y) @@ -500,8 +510,7 @@ def wrap_longitudes(lons): Longitudes wrapped into [-180:+180[ validity range """ - lons_wrap = (lons + 180) % (360) - 180 - return lons_wrap.astype(lons.dtype) + return (lons + 180) % 360 - 180 def recursive_dict_update(d, u): diff --git a/pyresample/version.py b/pyresample/version.py index 3e0309b..38c6ea6 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.7.0' +__version__ = '1.7.1' -- 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
