diff -Nru pyresample-1.5.0/appveyor.yml pyresample-1.7.0/appveyor.yml --- pyresample-1.5.0/appveyor.yml 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/appveyor.yml 2017-10-13 09:24:57.000000000 +0000 @@ -28,13 +28,13 @@ 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 @@ #on_success: # - TODO: upload the content of dist/*.whl to a public wheelhouse # - -notifications: - - provider: Slack - incoming_webhook: - secure: 7bEOYCIxHE5PkCF156zjxbJIeKkv7UpZulyn+Di09jqDlpvZjR0Qj3a1LK9AOjgwWgaQYbeI4BEYEdeq7pPxDs9snK8qvvbFDbuLAzg+HEQ= diff -Nru pyresample-1.5.0/.bumpversion.cfg pyresample-1.7.0/.bumpversion.cfg --- pyresample-1.5.0/.bumpversion.cfg 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/.bumpversion.cfg 2017-10-13 09:24:57.000000000 +0000 @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.5.0 +current_version = 1.7.0 commit = True tag = True diff -Nru pyresample-1.5.0/changelog.rst pyresample-1.7.0/changelog.rst --- pyresample-1.5.0/changelog.rst 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/changelog.rst 2017-10-13 09:24:57.000000000 +0000 @@ -2,6 +2,198 @@ ========= +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] +- Bump version: 1.6.0 ā†’ 1.6.1. [Martin Raspaud] +- Merge pull request #60 from pytroll/feature-dynamic-area. [David + Hoese] + + Add support for dynamic areas +- Merge branch 'develop' into feature-dynamic-area. [Martin Raspaud] +- Apply assert_allclose to proj dicts for tests. [Martin Raspaud] +- Fix some style issues. [Martin Raspaud] +- Set DynamicArea proj to `omerc` by default. [Martin Raspaud] +- Implement proposed changes in PR review. [Martin Raspaud] +- Use numpy's assert almost equal for area_extent comparisons. [Martin + Raspaud] +- Document the DynamicArea class. [Martin Raspaud] +- Fix optimal projection computation tests. [Martin Raspaud] +- Pep8 cleanup. [Martin Raspaud] +- Valid index computation optimization. [Martin Raspaud] +- Change bb computation api to use the whole proj_dict. [Martin Raspaud] +- Fix unittests for updated omerc computations. [Martin Raspaud] +- Use other azimuth direction for omerc. [Martin Raspaud] +- Flip x and y size in omerc projection. [Martin Raspaud] +- Bugfix typo. [Martin Raspaud] +- Allow lons and lats to be any array in bb computation. [Martin + Raspaud] +- Add SwathDefinition tests to the test suite. [Martin Raspaud] +- Support bounding box area computation from SwathDefintion. [Martin + Raspaud] + + This add support for computing a bounding box area from a swath definition that would fit optimally. The default projection is oblique mercator, with is optimal for locally received imager passes. +- Add support for dynamic areas. [Martin Raspaud] +- Merge pull request #70 from pytroll/feature-radius-parameters. [David + Hoese] + + Add 'proj4_radius_parameters' to calculate 'a' and 'b' from ellps +- Add tests for proj4_radius_parameters. [davidh-ssec] +- Fix typo in function call in radius parameters. [davidh-ssec] +- Add 'proj4_radius_parameters' to calculate 'a' and 'b' from ellps. + [davidh-ssec] +- Merge pull request #68 from pytroll/feature-56. [Martin Raspaud] + + Fix GridDefinition as permitted definition in preprocessing utils +- Add more preprocessing tests. [davidh-ssec] +- Fix preprocessing functions to use duck type on provided areas. + [davidh-ssec] +- Fix GridDefinition as permitted definition in preprocessing utils. + [davidh-ssec] + + +v1.6.0 (2017-09-12) +------------------- +- update changelog. [Martin Raspaud] +- Bump version: 1.5.0 ā†’ 1.6.0. [Martin Raspaud] +- Make sure x_size and y_size are ints. [Martin Raspaud] +- Merge pull request #69 from pytroll/bugfix-66. [Martin Raspaud] + + Fix write to mask affecting original mask in future versions of numpy + + Fixes #66 +- Add python 3.6 to travis tests. [davidh-ssec] +- Fix write to mask affecting original mask in future versions of numpy. + [davidh-ssec] + + Fix #66 + +- Merge pull request #67 from pytroll/bugfix-13. [Martin Raspaud] + + Rename `proj_x/y_coords` to `projection_x/y_coords` +- Rename `proj_x/y_coords` to `projection_x/y_coords` [davidh-ssec] + + Fix #13 + +- Merge pull request #63 from pytroll/feature-multiple-area-files. + [David Hoese] + + Parse multiple area files +- Fix tests_require in setup.py. [davidh-ssec] +- Use libgeos-dev to depend on the C++ libgeos-X.X.X and libgeos-c1. + [davidh-ssec] +- Add simple tests for parsing multiple yaml area strings. [davidh-ssec] +- Fix indentation in area file parsing functions. [davidh-ssec] +- Add ability to parse multiple area files at once. [davidh-ssec] +- Merge pull request #65 from pytroll/fix-numpy-1.13. [Martin Raspaud] + + Fix numpy 1.13 compatibility +- Fix boolean mask array usage in gaussian resampling. [davidh-ssec] + + In numpy 1.13 it is illegal to index an array with a boolean + array of a different size. + +- Add mock to test dependencies for python <3.3. [davidh-ssec] +- Use prepackaged numexpr in bdist_rpm. [Martin Raspaud] + + v1.5.0 (2017-05-02) ------------------- - update changelog. [Martin Raspaud] diff -Nru pyresample-1.5.0/debian/changelog pyresample-1.7.0/debian/changelog --- pyresample-1.5.0/debian/changelog 2017-09-05 08:22:41.000000000 +0000 +++ pyresample-1.7.0/debian/changelog 2017-10-14 07:46:10.000000000 +0000 @@ -1,14 +1,35 @@ -pyresample (1.5.0-3build2) artful; urgency=medium +pyresample (1.7.0-1) unstable; urgency=medium - * No-change rebuild to drop _PyFPE support. + * New upstream release + * Standard version bumped to v4.1.1 (no change) - -- Matthias Klose Tue, 05 Sep 2017 08:22:41 +0000 + -- Antonio Valentino Sat, 14 Oct 2017 07:46:10 +0000 -pyresample (1.5.0-3build1) artful; urgency=medium +pyresample (1.6.1-1) unstable; urgency=medium - * No-change rebuild to build to drop python3.5. + * New upstream release + * Standard version bumped to 4.1.0 (no change) + * Refresh all patches + * debian/control + - build-depend on python3-sphinx instead of python-sphinx + - remove un-necessary Testsuite: autopkgtest + * new debian/source.lintian-overrides file. + Suppress python-foo-but-no-python3-foo: the package only + contains common test data + + -- Antonio Valentino Sun, 24 Sep 2017 06:47:24 +0000 + +pyresample (1.6.0-1) unstable; urgency=medium + + * New upstream release + * debian/control + - drop dependenct from python3 mock (use the stdlib version) + - add dependency from six + * Refresh all patches + * debian/rules + - fix clean target - -- Matthias Klose Sat, 05 Aug 2017 18:39:16 +0000 + -- Antonio Valentino Sat, 16 Sep 2017 11:41:15 +0000 pyresample (1.5.0-3) unstable; urgency=medium diff -Nru pyresample-1.5.0/debian/control pyresample-1.7.0/debian/control --- pyresample-1.5.0/debian/control 2017-06-21 16:35:43.000000000 +0000 +++ pyresample-1.7.0/debian/control 2017-10-14 07:46:10.000000000 +0000 @@ -2,7 +2,6 @@ Maintainer: Debian GIS Project Uploaders: Antonio Valentino Section: python -Testsuite: autopkgtest Priority: optional Build-Depends: debhelper (>= 9.0.0), dh-python, @@ -18,7 +17,7 @@ python3-pyproj, python-configobj, python3-configobj, - python-sphinx, + python3-sphinx, python-numexpr, python3-numexpr, python-yaml, @@ -28,10 +27,11 @@ python-pil, python3-pil, python-mock, - python3-mock, + python-six, + python3-six, cython, cython3 -Standards-Version: 4.0.0 +Standards-Version: 4.1.1 Vcs-Browser: https://anonscm.debian.org/cgit/pkg-grass/pyresample.git Vcs-Git: https://anonscm.debian.org/git/pkg-grass/pyresample.git Homepage: https://github.com/pytroll/pyresample @@ -46,6 +46,7 @@ python-pyproj, python-yaml, python-pykdtree (>= 1.1.1), + python-six, python-pyresample-test Recommends: python-numexpr, python-pil, @@ -72,6 +73,7 @@ python3-pyproj, python3-yaml, python3-pykdtree (>= 1.1.1), + python3-six, python-pyresample-test Recommends: python3-numexpr, python3-pil, @@ -108,8 +110,7 @@ Package: python-pyresample-test Architecture: all Depends: ${misc:Depends}, - python-mock, - python3-mock + python-mock Description: Resampling of remote sensing data in Python (test suite) Pyresample is a Python package for resampling (reprojection) of earth observing satellite data. It handles both resampling of gridded data diff -Nru pyresample-1.5.0/debian/patches/0002-fix-proj4-initialization.patch pyresample-1.7.0/debian/patches/0002-fix-proj4-initialization.patch --- pyresample-1.5.0/debian/patches/0002-fix-proj4-initialization.patch 2017-06-04 16:06:33.000000000 +0000 +++ pyresample-1.7.0/debian/patches/0002-fix-proj4-initialization.patch 2017-09-23 18:11:13.000000000 +0000 @@ -21,10 +21,10 @@ YSIZE: 480 AREA_EXTENT: (-20037508.342789244, -10018754.171394622, 20037508.342789244, 10018754.171394622) diff --git a/pyresample/test/test_geometry.py b/pyresample/test/test_geometry.py -index 470ad52..5922c27 100644 +index 73a3050..98bd0a2 100644 --- a/pyresample/test/test_geometry.py +++ b/pyresample/test/test_geometry.py -@@ -408,7 +408,7 @@ class Test(unittest.TestCase): +@@ -309,7 +309,7 @@ class Test(unittest.TestCase): swath_def = geometry.SwathDefinition(lons, lats) filter_area = geometry.AreaDefinition('test', 'test', 'test', {'proj': 'eqc', 'lon_0': 0.0, @@ -33,7 +33,7 @@ 8, 8, (-20037508.34, -10018754.17, 20037508.34, 10018754.17)) filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0], -@@ -433,7 +433,7 @@ class Test(unittest.TestCase): +@@ -334,7 +334,7 @@ class Test(unittest.TestCase): data = np.array([1, 2, 3, 4]) filter_area = geometry.AreaDefinition('test', 'test', 'test', {'proj': 'eqc', 'lon_0': 0.0, @@ -42,7 +42,7 @@ 8, 8, (-20037508.34, -10018754.17, 20037508.34, 10018754.17)) filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0], -@@ -468,7 +468,7 @@ class Test(unittest.TestCase): +@@ -369,7 +369,7 @@ class Test(unittest.TestCase): data = np.dstack((data1, data2, data3)) filter_area = geometry.AreaDefinition('test', 'test', 'test', {'proj': 'eqc', 'lon_0': 0.0, @@ -51,7 +51,7 @@ 8, 8, (-20037508.34, -10018754.17, 20037508.34, 10018754.17)) filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0], -@@ -537,7 +537,7 @@ class Test(unittest.TestCase): +@@ -438,7 +438,7 @@ class Test(unittest.TestCase): def test_latlong_area(self): area_def = geometry.AreaDefinition('', '', '', diff -Nru pyresample-1.5.0/debian/rules pyresample-1.7.0/debian/rules --- pyresample-1.5.0/debian/rules 2017-06-25 19:36:51.000000000 +0000 +++ pyresample-1.7.0/debian/rules 2017-09-23 18:11:13.000000000 +0000 @@ -25,7 +25,7 @@ override_dh_auto_clean: dh_auto_clean --buildsystem=pybuild $(MAKE) -C docs clean - $(RM) -r build + $(RM) -r build pyresample.egg-info pyresample/ewa/*.so override_dh_installchangelogs: diff -Nru pyresample-1.5.0/debian/source.lintian-overrides pyresample-1.7.0/debian/source.lintian-overrides --- pyresample-1.5.0/debian/source.lintian-overrides 1970-01-01 00:00:00.000000000 +0000 +++ pyresample-1.7.0/debian/source.lintian-overrides 2017-09-24 17:35:48.000000000 +0000 @@ -0,0 +1,3 @@ +# the package only contains common test data +pyresample source: python-foo-but-no-python3-foo python-pyresample-test + diff -Nru pyresample-1.5.0/debian/tests/control pyresample-1.7.0/debian/tests/control --- pyresample-1.5.0/debian/tests/control 2017-06-04 16:06:33.000000000 +0000 +++ pyresample-1.7.0/debian/tests/control 2017-09-15 18:31:38.000000000 +0000 @@ -2,4 +2,4 @@ Depends: @builddeps@, python-pyresample, python-pyresample-test, python-all, python-mock Tests: python3 -Depends: @builddeps@, python3-pyresample, python-pyresample-test, python3-all, python3-mock +Depends: @builddeps@, python3-pyresample, python-pyresample-test, python3-all diff -Nru pyresample-1.5.0/docs/source/geo_def.rst pyresample-1.7.0/docs/source/geo_def.rst --- pyresample-1.5.0/docs/source/geo_def.rst 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/docs/source/geo_def.rst 2017-10-13 09:24:57.000000000 +0000 @@ -267,13 +267,18 @@ ------------------------------------- A ***definition** object allows for retrieval of geographic coordinates using array slicing (slice stepping is currently not supported). -All ***definition** objects exposes the coordinates **lons**, **lats** and **cartesian_coords**. +All ***definition** objects expose the coordinates **lons**, **lats** and **cartesian_coords**. AreaDefinition exposes the full set of projection coordinates as **projection_x_coords** and **projection_y_coords**. Note that in the case of projection coordinates expressed in longitude and latitude, **projection_x_coords** will be longitude and **projection_y_coords** will be latitude. +.. versionchanged:: 1.5.1 + + Renamed `proj_x_coords` to `projection_x_coords` and `proj_y_coords` + to `projection_y_coords`. + Get full coordinate set: .. doctest:: @@ -306,8 +311,8 @@ ... x_size, y_size, area_extent) >>> cart_subset = area_def.get_cartesian_coords()[100:200, 350:] -If only the 1D range of a projection coordinate is required it can be extraxted -using the **proj_x_coord** or **proj_y_coords** property of a geographic coordinate +If only the 1D range of a projection coordinate is required it can be extracted +using the **projection_x_coord** or **projection_y_coords** property of a geographic coordinate .. doctest:: @@ -321,7 +326,7 @@ >>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625) >>> area_def = utils.get_area_def(area_id, description, proj_id, projection, ... x_size, y_size, area_extent) - >>> proj_x_range = area_def.proj_x_coords + >>> proj_x_range = area_def.projection_x_coords Spherical geometry operations ----------------------------- diff -Nru pyresample-1.5.0/docs/source/swath.rst pyresample-1.7.0/docs/source/swath.rst --- pyresample-1.5.0/docs/source/swath.rst 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/docs/source/swath.rst 2017-10-13 09:24:57.000000000 +0000 @@ -8,7 +8,7 @@ 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 @@ >>> 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 @@ ... 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 @@ * **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 @@ **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 @@ >>> 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 -Nru pyresample-1.5.0/pyresample/bilinear/__init__.py pyresample-1.7.0/pyresample/bilinear/__init__.py --- pyresample-1.5.0/pyresample/bilinear/__init__.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/bilinear/__init__.py 2017-10-13 09:24:57.000000000 +0000 @@ -30,6 +30,7 @@ import numpy as np from pyproj import Proj +import warnings from pyresample import kd_tree @@ -97,6 +98,12 @@ 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 @@ # 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 @@ 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 -Nru pyresample-1.5.0/pyresample/geometry.py pyresample-1.7.0/pyresample/geometry.py --- pyresample-1.5.0/pyresample/geometry.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/geometry.py 2017-10-13 09:24:57.000000000 +0000 @@ -24,12 +24,16 @@ import warnings from collections import OrderedDict +from logging import getLogger import numpy as np import yaml +from pyproj import Geod, Proj from pyresample import _spatial_mp, utils +logger = getLogger(__name__) + class DimensionError(Exception): pass @@ -59,6 +63,8 @@ if type(lons) != type(lats): raise TypeError('lons and lats must be of same type') elif lons is not None: + lons = np.asanyarray(lons) + lats = np.asanyarray(lats) if lons.shape != lats.shape: raise ValueError('lons and lats must have same shape') @@ -351,6 +357,8 @@ """Base class for geometry definitions defined by lons and lats only""" def __init__(self, lons, lats, nprocs=1): + lons = np.asanyarray(lons) + lats = np.asanyarray(lats) super(CoordinateDefinition, self).__init__(lons, lats, nprocs) if lons.shape == lats.shape and lons.dtype == lats.dtype: self.shape = lons.shape @@ -449,12 +457,183 @@ """ def __init__(self, lons, lats, nprocs=1): + lons = np.asanyarray(lons) + lats = np.asanyarray(lats) super(SwathDefinition, self).__init__(lons, lats, nprocs) if lons.shape != lats.shape: raise ValueError('lon and lat arrays must have same shape') elif lons.ndim > 2: raise ValueError('Only 1 and 2 dimensional swaths are allowed') + def 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 + lon1, lon2 = np.asanyarray(self.lons[[0, -1], int(cols / 2)]) + lat1, lat, lat2 = np.asanyarray( + self.lats[[0, int(lines / 2), -1], int(cols / 2)]) + + proj_dict2points = {'proj': 'omerc', 'lat_0': lat, 'ellps': ellipsoid, + 'lat_1': lat1, 'lon_1': lon1, + 'lat_2': lat2, 'lon_2': lon2} + + lonc, lat0 = Proj(**proj_dict2points)(0, 0, inverse=True) + az1, az2, dist = Geod(**proj_dict2points).inv(lonc, lat0, lon1, lat1) + del az2, dist + return {'proj': 'omerc', 'alpha': float(az1), + 'lat_0': float(lat0), 'lonc': float(lonc), + 'no_rot': True, 'ellps': ellipsoid} + + def get_edge_lonlats(self): + """Get the concatenated boundary of the current swath.""" + lons, lats = self.get_boundary_lonlats() + blons = np.ma.concatenate([lons.side1, lons.side2, + lons.side3, lons.side4]) + blats = np.ma.concatenate([lats.side1, lats.side2, + lats.side3, lats.side4]) + return blons, blats + + def compute_bb_proj_params(self, proj_dict): + projection = proj_dict['proj'] + ellipsoid = proj_dict.get('ellps', 'WGS84') + if projection == 'omerc': + return self._compute_omerc_parameters(ellipsoid) + else: + raise NotImplementedError('Only omerc supported for now.') + + def compute_optimal_bb_area(self, proj_dict=None): + """Compute the "best" bounding box area for this swath with `proj_dict`. + + By default, the projection is Oblique Mercator (`omerc` in proj.4), in + which case the right projection angle `alpha` is computed from the + swath centerline. For other projections, only the appropriate center of + projection and area extents are computed. + """ + if proj_dict is None: + proj_dict = {} + projection = proj_dict.setdefault('proj', 'omerc') + area_id = projection + '_otf' + description = 'On-the-fly ' + projection + ' area' + lines, cols = self.lons.shape + x_size = int(cols * 1.1) + y_size = int(lines * 1.1) + + proj_dict = self.compute_bb_proj_params(proj_dict) + + if projection == 'omerc': + x_size, y_size = y_size, x_size + + area = DynamicAreaDefinition(area_id, description, proj_dict) + lons, lats = self.get_edge_lonlats() + return area.freeze((lons, lats), size=(x_size, y_size)) + + +class DynamicAreaDefinition(object): + """An AreaDefintion containing just a subset of the needed parameters. + + The purpose of this class is to be able to adapt the area extent and size of + the area to a given set of longitudes and latitudes, such that e.g. polar + satellite granules can be resampled optimaly to a give projection.""" + + def __init__(self, area_id=None, description=None, proj_dict=None, + x_size=None, y_size=None, area_extent=None, + optimize_projection=False): + """Initialize the DynamicAreaDefinition. + + area_id: + The name of the area. + description: + The description of the area. + proj_dict: + The dictionary of projection parameters. Doesn't have to be complete. + x_size, y_size: + The size of the resulting area. + area_extent: + The area extent of the area. + optimize_projection: + Whether the projection parameters have to be optimized. + """ + self.area_id = area_id + self.description = description + self.proj_dict = proj_dict + self.x_size = x_size + self.y_size = y_size + self.area_extent = area_extent + self.optimize_projection = optimize_projection + + def compute_domain(self, corners, resolution=None, size=None): + """Compute size and area_extent from corners and [size or resolution] + info.""" + if resolution is not None and size is not None: + raise ValueError("Both resolution and size can't be provided.") + + if size: + x_size, y_size = size + x_resolution = (corners[2] - corners[0]) * 1.0 / (x_size - 1) + y_resolution = (corners[3] - corners[1]) * 1.0 / (y_size - 1) + + if resolution: + try: + x_resolution, y_resolution = resolution + except TypeError: + x_resolution = y_resolution = resolution + x_size = int(np.rint((corners[2] - corners[0]) * 1.0 / + x_resolution + 1)) + y_size = int(np.rint((corners[3] - corners[1]) * 1.0 / + y_resolution + 1)) + + area_extent = (corners[0] - x_resolution / 2, + corners[1] - y_resolution / 2, + corners[2] + x_resolution / 2, + corners[3] + y_resolution / 2) + return area_extent, x_size, y_size + + def freeze(self, lonslats=None, + resolution=None, size=None, + proj_info=None): + """Create an AreaDefintion from this area with help of some extra info. + + lonlats: + the geographical coordinates to contain in the resulting area. + resolution: + the resolution of the resulting area. + size: + the size of the resulting area. + proj_info: + complementing parameters to the projection info. + + Resolution and Size parameters are ignored if the instance is created + with the `optimize_projection` flag set to True. + """ + if proj_info is not None: + self.proj_dict.update(proj_info) + + if self.optimize_projection: + return lonslats.compute_optimal_bb_area(self.proj_dict) + + if not self.area_extent or not self.x_size or not self.y_size: + proj4 = Proj(**self.proj_dict) + try: + lons, lats = lonslats + except (TypeError, ValueError): + lons, lats = lonslats.get_lonlats() + xarr, yarr = proj4(np.asarray(lons), np.asarray(lats)) + corners = [np.min(xarr), np.min(yarr), np.max(xarr), np.max(yarr)] + + domain = self.compute_domain(corners, resolution, size) + self.area_extent, self.x_size, self.y_size = domain + + return AreaDefinition(self.area_id, self.description, '', + self.proj_dict, self.x_size, self.y_size, + self.area_extent) + class AreaDefinition(BaseDefinition): @@ -541,8 +720,8 @@ self.area_id = area_id self.name = name self.proj_id = proj_id - self.x_size = x_size - self.y_size = y_size + self.x_size = int(x_size) + self.y_size = int(y_size) self.shape = (y_size, x_size) if lons is not None: if lons.shape != self.shape: @@ -575,11 +754,15 @@ self.pixel_offset_x = -self.area_extent[0] / self.pixel_size_x self.pixel_offset_y = self.area_extent[3] / self.pixel_size_y - self.projection_x_coords = None - self.projection_y_coords = None + self._projection_x_coords = None + self._projection_y_coords = None 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 @@ -588,7 +771,7 @@ 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) @@ -659,8 +842,8 @@ (see get_lonlats). """ p = _spatial_mp.Proj(self.proj4_string) - x = self.proj_x_coords - y = self.proj_y_coords + x = self.projection_x_coords + y = self.projection_y_coords return p(y[y.size - cols], x[x.size - rows], inverse=True) def lonlat2colrow(self, lons, lats): @@ -780,6 +963,34 @@ 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 @@ -788,7 +999,7 @@ 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 ------- @@ -807,12 +1018,12 @@ else: return val - if self.projection_x_coords is not None and self.projection_y_coords is not None: + if self._projection_x_coords is not None and self._projection_y_coords is not None: # Projection coords are cached if data_slice is None: - return self.projection_x_coords, self.projection_y_coords + return self._projection_x_coords, self._projection_y_coords else: - return self.projection_x_coords[data_slice], self.projection_y_coords[data_slice] + return self._projection_x_coords[data_slice], self._projection_y_coords[data_slice] is_single_value = False is_1d_select = False @@ -870,17 +1081,9 @@ 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 @@ -893,20 +1096,32 @@ if cache and data_slice is None: # Cache the result if requested - self.projection_x_coords = target_x - self.projection_y_coords = target_y + self._projection_x_coords = target_x + self._projection_y_coords = target_y return target_x, target_y @property - def proj_x_coords(self): + def projection_x_coords(self): return self.get_proj_coords(data_slice=(0, slice(None)))[0] @property - def proj_y_coords(self): + def projection_y_coords(self): return self.get_proj_coords(data_slice=(slice(None), 0))[1] @property + def proj_x_coords(self): + 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) + return self.projection_y_coords + + @property def outer_boundary_corners(self): """Returns the lon,lat of the outer edges of the corner points """ @@ -923,6 +1138,39 @@ 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. @@ -1026,8 +1274,8 @@ if axis == 0: area_extent = combine_area_extents_vertical(area1, area2) - x_size = area1.x_size - y_size = area1.y_size + area2.y_size + x_size = int(area1.x_size) + y_size = int(area1.y_size + area2.y_size) else: raise NotImplementedError('Only vertical contatenation is supported.') return AreaDefinition(area1.area_id, area1.name, area1.proj_id, @@ -1112,6 +1360,20 @@ 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 -Nru pyresample-1.5.0/pyresample/image.py pyresample-1.7.0/pyresample/image.py --- pyresample-1.5.0/pyresample/image.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/image.py 2017-10-13 09:24:57.000000000 +0000 @@ -21,32 +21,32 @@ 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 @@ ---------- 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 @@ ---------- 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 @@ ---------- 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 @@ 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 @@ 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 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 @@ 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 @@ 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 @@ 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 @@ 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 -Nru pyresample-1.5.0/pyresample/kd_tree.py pyresample-1.7.0/pyresample/kd_tree.py --- pyresample-1.5.0/pyresample/kd_tree.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/kd_tree.py 2017-10-13 09:24:57.000000000 +0000 @@ -21,13 +21,22 @@ from __future__ import absolute_import +import sys import types import warnings -import sys +from logging import getLogger import numpy as np -from pyresample import geometry, data_reduce, _spatial_mp +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 @@ -66,20 +75,20 @@ ---------- source_geo_def : object Geometry definition of source - data : numpy array + data : numpy array 1d array of single channel data points or (source_size, k) array of k channels of datapoints target_geo_def : object Geometry definition of target - radius_of_influence : float + radius_of_influence : float Cut off distance in meters epsilon : float, optional Allowed uncertainty in meters. Increasing uncertainty reduces execution time fill_value : int or None, optional Set undetermined pixels to this value. - If fill_value is None a masked array is returned - with undetermined pixels masked + If fill_value is None a masked array is returned + with undetermined pixels masked reduce_data : bool, optional Perform initial coarse reduction of source dataset in order to reduce execution time @@ -91,7 +100,7 @@ Returns ------- - data : numpy array + data : numpy array Source data resampled to target geometry """ @@ -110,26 +119,26 @@ ---------- source_geo_def : object Geometry definition of source - data : numpy array + data : numpy array Array of single channel data points or (source_geo_def.shape, k) array of k channels of datapoints target_geo_def : object Geometry definition of target - radius_of_influence : float + radius_of_influence : float Cut off distance in meters - sigmas : list of floats or float - List of sigmas to use for the gauss weighting of each + sigmas : list of floats or float + List of sigmas to use for the gauss weighting of each channel 1 to k, w_k = exp(-dist^2/sigma_k^2). If only one channel is resampled sigmas is a single float value. - neighbours : int, optional + neighbours : int, optional The number of neigbours to consider for each grid point epsilon : float, optional Allowed uncertainty in meters. Increasing uncertainty reduces execution time - fill_value : {int, None}, optional + fill_value : {int, None}, optional Set undetermined pixels to this value. - If fill_value is None a masked array is returned - with undetermined pixels masked + If fill_value is None a masked array is returned + with undetermined pixels masked reduce_data : bool, optional Perform initial coarse reduction of source dataset in order to reduce execution time @@ -148,7 +157,7 @@ data, stddev, counts : numpy array, numpy array, numpy array (if with_uncert == True) Source data resampled to target geometry. Weighted standard devaition for all pixels having more than one source value - Counts of number of source values used in weighting per pixel + Counts of number of source values used in weighting per pixel """ def gauss(sigma): @@ -190,27 +199,27 @@ ---------- source_geo_def : object Geometry definition of source - data : numpy array + data : numpy array Array of single channel data points or (source_geo_def.shape, k) array of k channels of datapoints target_geo_def : object Geometry definition of target - radius_of_influence : float + radius_of_influence : float Cut off distance in meters - weight_funcs : list of function objects or function object - List of weight functions f(dist) to use for the weighting + weight_funcs : list of function objects or function object + List of weight functions f(dist) to use for the weighting of each channel 1 to k. If only one channel is resampled weight_funcs is a single function object. - neighbours : int, optional + neighbours : int, optional The number of neigbours to consider for each grid point epsilon : float, optional Allowed uncertainty in meters. Increasing uncertainty reduces execution time - fill_value : {int, None}, optional + fill_value : {int, None}, optional Set undetermined pixels to this value. - If fill_value is None a masked array is returned - with undetermined pixels masked + If fill_value is None a masked array is returned + with undetermined pixels masked reduce_data : bool, optional Perform initial coarse reduction of source dataset in order to reduce execution time @@ -282,9 +291,9 @@ Geometry definition of source target_geo_def : object Geometry definition of target - radius_of_influence : float + radius_of_influence : float Cut off distance in meters - neighbours : int, optional + neighbours : int, optional The number of neigbours to consider for each grid point epsilon : float, optional Allowed uncertainty in meters. Increasing uncertainty @@ -300,7 +309,7 @@ Returns ------- - (valid_input_index, valid_output_index, + (valid_input_index, valid_output_index, index_array, distance_array) : tuple of numpy arrays Neighbour resampling info """ @@ -390,8 +399,8 @@ """Find indices of reduced inputput data""" source_lons, source_lats = source_geo_def.get_lonlats(nprocs=nprocs) - source_lons = source_lons.ravel() - source_lats = source_lats.ravel() + source_lons = np.asanyarray(source_lons).ravel() + source_lats = np.asanyarray(source_lats).ravel() if source_lons.size == 0 or source_lats.size == 0: raise ValueError('Cannot resample empty data set') @@ -400,9 +409,8 @@ raise ValueError('Mismatch between lons and lats') # Remove illegal values - valid_data = ((source_lons >= -180) & (source_lons <= 180) & - (source_lats <= 90) & (source_lats >= -90)) - valid_input_index = np.ones(source_geo_def.size, dtype=np.bool) + valid_input_index = ((source_lons >= -180) & (source_lons <= 180) & + (source_lats <= 90) & (source_lats >= -90)) if reduce_data: # Reduce dataset @@ -415,16 +423,15 @@ geometry.AreaDefinition))): # Resampling from swath to grid or from grid to grid lonlat_boundary = target_geo_def.get_boundary_lonlats() - valid_input_index = \ + + # Combine reduced and legal values + valid_input_index &= \ data_reduce.get_valid_index_from_lonlat_boundaries( lonlat_boundary[0], lonlat_boundary[1], source_lons, source_lats, radius_of_influence) - # Combine reduced and legal values - valid_input_index = (valid_data & valid_input_index) - if(isinstance(valid_input_index, np.ma.core.MaskedArray)): # Make sure valid_input_index is not a masked array valid_input_index = valid_input_index.filled(False) @@ -469,7 +476,7 @@ """ if not isinstance(source_geo_def, geometry.BaseDefinition): raise TypeError('source_geo_def must be of geometry type') - + #Get reduced cartesian coordinates and flatten them source_cartesian_coords = source_geo_def.get_cartesian_coords(nprocs=nprocs) input_coords = geometry._flatten_cartesian_coords(source_cartesian_coords) @@ -599,20 +606,20 @@ distance_array : numpy array, optional distance_array from get_neighbour_info Not needed for 'nn' resample type - weight_funcs : list of function objects or function object, optional - List of weight functions f(dist) to use for the weighting + weight_funcs : list of function objects or function object, optional + List of weight functions f(dist) to use for the weighting of each channel 1 to k. If only one channel is resampled weight_funcs is a single function object. Must be supplied when using 'custom' resample type fill_value : int or None, optional Set undetermined pixels to this value. - If fill_value is None a masked array is returned + If fill_value is None a masked array is returned with undetermined pixels masked Returns ------- - result : numpy array + result : numpy array Source data resampled to target geometry """ @@ -716,7 +723,7 @@ # Get nearest neighbour using array indexing index_mask = (index_array == input_size) new_index_array = np.where(index_mask, 0, index_array) - result = new_data[new_index_array] + result = new_data[new_index_array].copy() result[index_mask] = fill_value else: # Calculate result using weighting. @@ -807,11 +814,22 @@ # Calculate final stddev new_valid_index = (count > 1) - v1 = norm[new_valid_index] - v2 = norm_sqr[new_valid_index] - stddev[new_valid_index] = np.sqrt( - (v1 / (v1 ** 2 - v2)) * stddev[new_valid_index]) - stddev[~new_valid_index] = np.NaN + if stddev.ndim >= 2: + # If given more than 1 input data array + new_valid_index = new_valid_index[:, 0] + for i in range(stddev.shape[-1]): + v1 = norm[new_valid_index, i] + v2 = norm_sqr[new_valid_index, i] + stddev[new_valid_index, i] = np.sqrt( + (v1 / (v1 ** 2 - v2)) * stddev[new_valid_index, i]) + stddev[~new_valid_index, i] = np.NaN + else: + # If given single input data array + v1 = norm[new_valid_index] + v2 = norm_sqr[new_valid_index] + stddev[new_valid_index] = np.sqrt( + (v1 / (v1 ** 2 - v2)) * stddev[new_valid_index]) + stddev[~new_valid_index] = np.NaN # Add fill values result[np.invert(result_valid_index)] = fill_value @@ -868,6 +886,274 @@ 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 -Nru pyresample-1.5.0/pyresample/test/test_bilinear.py pyresample-1.7.0/pyresample/test/test_bilinear.py --- pyresample-1.5.0/pyresample/test/test_bilinear.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/test_bilinear.py 2017-10-13 09:24:57.000000000 +0000 @@ -207,7 +207,7 @@ 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 @@ 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 -Nru pyresample-1.5.0/pyresample/test/test_ewa_ll2cr.py pyresample-1.7.0/pyresample/test/test_ewa_ll2cr.py --- pyresample-1.5.0/pyresample/test/test_ewa_ll2cr.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/test_ewa_ll2cr.py 2017-10-13 09:24:57.000000000 +0000 @@ -24,6 +24,7 @@ import sys import logging import numpy as np +from pyresample.test.utils import create_test_longitude, create_test_latitude if sys.version_info < (2, 7): import unittest2 as unittest else: @@ -32,28 +33,6 @@ LOG = logging.getLogger(__name__) -def create_test_longitude(start, stop, shape, twist_factor=0.0, dtype=np.float32): - if start > stop: - stop += 360.0 - - lon_row = np.linspace(start, stop, num=shape[1]).astype(dtype) - twist_array = np.arange(shape[0]).reshape((shape[0], 1)) * twist_factor - lon_array = np.repeat([lon_row], shape[0], axis=0) - lon_array += twist_array - - if stop > 360.0: - lon_array[lon_array > 360.0] -= 360 - return lon_array - - -def create_test_latitude(start, stop, shape, twist_factor=0.0, dtype=np.float32): - lat_col = np.linspace(start, stop, num=shape[0]).astype(dtype).reshape((shape[0], 1)) - twist_array = np.arange(shape[1]) * twist_factor - lat_array = np.repeat(lat_col, shape[1], axis=1) - lat_array += twist_array - return lat_array - - dynamic_wgs84 = { "grid_name": "test_wgs84_fit", "origin_x": None, diff -Nru pyresample-1.5.0/pyresample/test/test_geometry.py pyresample-1.7.0/pyresample/test/test_geometry.py --- pyresample-1.5.0/pyresample/test/test_geometry.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/test_geometry.py 2017-10-13 09:24:57.000000000 +0000 @@ -4,7 +4,6 @@ import sys import numpy as np -from mock import MagicMock, patch from pyresample import geo_filter, geometry from pyresample.geometry import (IncompatibleAreas, @@ -12,6 +11,12 @@ concatenate_area_defs) from pyresample.test.utils import catch_warnings +try: + from unittest.mock import MagicMock, patch +except ImportError: + # separate mock package py<3.3 + from mock import MagicMock, patch + if sys.version_info < (2, 7): import unittest2 as unittest else: @@ -160,39 +165,6 @@ "BaseDefinition did not maintain dtype of longitudes (in:%s out:%s)" % (lons2_ints.dtype, lons.dtype,)) - def test_swath(self): - lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100)) - lats1 = np.fromfunction( - lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100)) - - swath_def = geometry.SwathDefinition(lons1, lats1) - - lons2, lats2 = swath_def.get_lonlats() - - self.assertFalse(id(lons1) != id(lons2) or id(lats1) != id(lats2), - msg='Caching of swath coordinates failed') - - def test_swath_wrap(self): - lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100)) - lats1 = np.fromfunction( - lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100)) - - lons1 += 180. - with catch_warnings() as w1: - swath_def = geometry.BaseDefinition(lons1, lats1) - self.assertFalse( - len(w1) != 1, 'Failed to trigger a warning on longitude wrapping') - self.assertFalse(('-180:+180' not in str(w1[0].message)), - 'Failed to trigger correct warning about longitude wrapping') - - lons2, lats2 = swath_def.get_lonlats() - - self.assertTrue(id(lons1) != id(lons2), - msg='Caching of swath coordinates failed with longitude wrapping') - - self.assertTrue(lons2.min() > -180 and lons2.max() < 180, - 'Wrapping of longitudes failed for SwathDefinition') - def test_area_equal(self): area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', {'a': '6378144.0', @@ -255,23 +227,6 @@ self.assertFalse( area_def == msg_area, 'area_defs are not expected to be equal') - def test_swath_equal(self): - lons = np.array([1.2, 1.3, 1.4, 1.5]) - lats = np.array([65.9, 65.86, 65.82, 65.78]) - swath_def = geometry.SwathDefinition(lons, lats) - swath_def2 = geometry.SwathDefinition(lons, lats) - self.assertFalse( - swath_def != swath_def2, 'swath_defs are not equal as expected') - - def test_swath_not_equal(self): - lats1 = np.array([65.9, 65.86, 65.82, 65.78]) - lons = np.array([1.2, 1.3, 1.4, 1.5]) - lats2 = np.array([65.91, 65.85, 65.80, 65.75]) - swath_def = geometry.SwathDefinition(lons, lats1) - swath_def2 = geometry.SwathDefinition(lons, lats2) - self.assertFalse( - swath_def == swath_def2, 'swath_defs are not expected to be equal') - def test_swath_equal_area(self): area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD', {'a': '6378144.0', @@ -348,60 +303,6 @@ self.assertFalse( area_def == swath_def, "swath_def and area_def should be different") - def test_concat_1d(self): - lons1 = np.array([1, 2, 3]) - lats1 = np.array([1, 2, 3]) - lons2 = np.array([4, 5, 6]) - lats2 = np.array([4, 5, 6]) - swath_def1 = geometry.SwathDefinition(lons1, lats1) - swath_def2 = geometry.SwathDefinition(lons2, lats2) - swath_def_concat = swath_def1.concatenate(swath_def2) - expected = np.array([1, 2, 3, 4, 5, 6]) - self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and - np.array_equal(swath_def_concat.lons, expected), - 'Failed to concatenate 1D swaths') - - def test_concat_2d(self): - lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) - lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) - lons2 = np.array([[4, 5, 6], [6, 7, 8]]) - lats2 = np.array([[4, 5, 6], [6, 7, 8]]) - swath_def1 = geometry.SwathDefinition(lons1, lats1) - swath_def2 = geometry.SwathDefinition(lons2, lats2) - swath_def_concat = swath_def1.concatenate(swath_def2) - expected = np.array( - [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]]) - self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and - np.array_equal(swath_def_concat.lons, expected), - 'Failed to concatenate 2D swaths') - - def test_append_1d(self): - lons1 = np.array([1, 2, 3]) - lats1 = np.array([1, 2, 3]) - lons2 = np.array([4, 5, 6]) - lats2 = np.array([4, 5, 6]) - swath_def1 = geometry.SwathDefinition(lons1, lats1) - swath_def2 = geometry.SwathDefinition(lons2, lats2) - swath_def1.append(swath_def2) - expected = np.array([1, 2, 3, 4, 5, 6]) - self.assertTrue(np.array_equal(swath_def1.lons, expected) and - np.array_equal(swath_def1.lons, expected), - 'Failed to append 1D swaths') - - def test_append_2d(self): - lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) - lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) - lons2 = np.array([[4, 5, 6], [6, 7, 8]]) - lats2 = np.array([[4, 5, 6], [6, 7, 8]]) - swath_def1 = geometry.SwathDefinition(lons1, lats1) - swath_def2 = geometry.SwathDefinition(lons2, lats2) - swath_def1.append(swath_def2) - expected = np.array( - [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]]) - self.assertTrue(np.array_equal(swath_def1.lons, expected) and - np.array_equal(swath_def1.lons, expected), - 'Failed to append 2D swaths') - def test_grid_filter_valid(self): lons = np.array([-170, -30, 30, 170]) lats = np.array([20, -40, 50, -80]) @@ -505,7 +406,7 @@ -909968.64000000001, 1029087.28, 1490031.3600000001]) - proj_x_boundary, proj_y_boundary = area_def.proj_x_coords, area_def.proj_y_coords + proj_x_boundary, proj_y_boundary = area_def.projection_x_coords, area_def.projection_y_coords expected_x = np.array([-1250912.72, -1010912.72, -770912.72, -530912.72, -290912.72, -50912.72, 189087.28, 429087.28, 669087.28, 909087.28]) @@ -701,6 +602,205 @@ self.assertTrue((y__.data == y_expects).all()) +def assert_np_dict_allclose(dict1, dict2): + + assert set(dict1.keys()) == set(dict2.keys()) + for key, val in dict1.items(): + try: + np.testing.assert_allclose(val, dict2[key]) + except TypeError: + assert(val == dict2[key]) + + +class TestSwathDefinition(unittest.TestCase): + """Test the SwathDefinition.""" + + def test_swath(self): + lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100)) + lats1 = np.fromfunction( + lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100)) + + swath_def = geometry.SwathDefinition(lons1, lats1) + + lons2, lats2 = swath_def.get_lonlats() + + self.assertFalse(id(lons1) != id(lons2) or id(lats1) != id(lats2), + msg='Caching of swath coordinates failed') + + def test_swath_wrap(self): + lons1 = np.fromfunction(lambda y, x: 3 + (10.0 / 100) * x, (5000, 100)) + lats1 = np.fromfunction( + lambda y, x: 75 - (50.0 / 5000) * y, (5000, 100)) + + lons1 += 180. + with catch_warnings() as w1: + swath_def = geometry.BaseDefinition(lons1, lats1) + self.assertFalse( + len(w1) != 1, 'Failed to trigger a warning on longitude wrapping') + self.assertFalse(('-180:+180' not in str(w1[0].message)), + 'Failed to trigger correct warning about longitude wrapping') + + lons2, lats2 = swath_def.get_lonlats() + + self.assertTrue(id(lons1) != id(lons2), + msg='Caching of swath coordinates failed with longitude wrapping') + + self.assertTrue(lons2.min() > -180 and lons2.max() < 180, + 'Wrapping of longitudes failed for SwathDefinition') + + def test_concat_1d(self): + lons1 = np.array([1, 2, 3]) + lats1 = np.array([1, 2, 3]) + lons2 = np.array([4, 5, 6]) + lats2 = np.array([4, 5, 6]) + swath_def1 = geometry.SwathDefinition(lons1, lats1) + swath_def2 = geometry.SwathDefinition(lons2, lats2) + swath_def_concat = swath_def1.concatenate(swath_def2) + expected = np.array([1, 2, 3, 4, 5, 6]) + self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and + np.array_equal(swath_def_concat.lons, expected), + 'Failed to concatenate 1D swaths') + + def test_concat_2d(self): + lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) + lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) + lons2 = np.array([[4, 5, 6], [6, 7, 8]]) + lats2 = np.array([[4, 5, 6], [6, 7, 8]]) + swath_def1 = geometry.SwathDefinition(lons1, lats1) + swath_def2 = geometry.SwathDefinition(lons2, lats2) + swath_def_concat = swath_def1.concatenate(swath_def2) + expected = np.array( + [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]]) + self.assertTrue(np.array_equal(swath_def_concat.lons, expected) and + np.array_equal(swath_def_concat.lons, expected), + 'Failed to concatenate 2D swaths') + + def test_append_1d(self): + lons1 = np.array([1, 2, 3]) + lats1 = np.array([1, 2, 3]) + lons2 = np.array([4, 5, 6]) + lats2 = np.array([4, 5, 6]) + swath_def1 = geometry.SwathDefinition(lons1, lats1) + swath_def2 = geometry.SwathDefinition(lons2, lats2) + swath_def1.append(swath_def2) + expected = np.array([1, 2, 3, 4, 5, 6]) + self.assertTrue(np.array_equal(swath_def1.lons, expected) and + np.array_equal(swath_def1.lons, expected), + 'Failed to append 1D swaths') + + def test_append_2d(self): + lons1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) + lats1 = np.array([[1, 2, 3], [3, 4, 5], [5, 6, 7]]) + lons2 = np.array([[4, 5, 6], [6, 7, 8]]) + lats2 = np.array([[4, 5, 6], [6, 7, 8]]) + swath_def1 = geometry.SwathDefinition(lons1, lats1) + swath_def2 = geometry.SwathDefinition(lons2, lats2) + swath_def1.append(swath_def2) + expected = np.array( + [[1, 2, 3], [3, 4, 5], [5, 6, 7], [4, 5, 6], [6, 7, 8]]) + self.assertTrue(np.array_equal(swath_def1.lons, expected) and + np.array_equal(swath_def1.lons, expected), + 'Failed to append 2D swaths') + + def test_swath_equal(self): + """Test swath equality.""" + lons = np.array([1.2, 1.3, 1.4, 1.5]) + lats = np.array([65.9, 65.86, 65.82, 65.78]) + swath_def = geometry.SwathDefinition(lons, lats) + swath_def2 = geometry.SwathDefinition(lons, lats) + self.assertFalse( + swath_def != swath_def2, 'swath_defs are not equal as expected') + + def test_swath_not_equal(self): + """Test swath inequality.""" + lats1 = np.array([65.9, 65.86, 65.82, 65.78]) + lons = np.array([1.2, 1.3, 1.4, 1.5]) + lats2 = np.array([65.91, 65.85, 65.80, 65.75]) + swath_def = geometry.SwathDefinition(lons, lats1) + swath_def2 = geometry.SwathDefinition(lons, lats2) + self.assertFalse( + swath_def == swath_def2, 'swath_defs are not expected to be equal') + + def test_compute_omerc_params(self): + """Test omerc parameters computation.""" + lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469], + [80.84000396728516, 60.74200439453125, 34.08500289916992], + [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T + + lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906], + [79.11000061035156, 7.284000396728516, -5.107000350952148], + [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T + + area = geometry.SwathDefinition(lons, lats) + proj_dict = {'no_rot': True, 'lonc': 5.340645620216994, + 'ellps': 'WGS84', 'proj': 'omerc', + 'alpha': 19.022450179020247, 'lat_0': 60.7420043944989} + assert_np_dict_allclose(area._compute_omerc_parameters('WGS84'), + proj_dict) + + def test_get_edge_lonlats(self): + """Test the `get_edge_lonlats` functionality.""" + lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469], + [80.84000396728516, 60.74200439453125, 34.08500289916992], + [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T + + lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906], + [79.11000061035156, 7.284000396728516, -5.107000350952148], + [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T + + area = geometry.SwathDefinition(lons, lats) + lons, lats = area.get_edge_lonlats() + + np.testing.assert_allclose(lons, [-90.67900085, 79.11000061, 81.26400757, + 81.26400757, 29.67200089, 10.26000023, + 10.26000023, -5.10700035, -21.52500153, + -21.52500153, -21.56500053, -90.67900085]) + np.testing.assert_allclose(lats, [85.23900604, 80.84000397, 67.07600403, + 67.07600403, 54.14700317, 30.54700089, + 30.54700089, 34.0850029, 35.58000183, + 35.58000183, 62.25600433, 85.23900604]) + + lats = np.array([[80., 80., 80.], + [80., 90., 80], + [80., 80., 80.]]).T + + lons = np.array([[-45., 0., 45.], + [-90, 0., 90.], + [-135., 180., 135.]]).T + + area = geometry.SwathDefinition(lons, lats) + lons, lats = area.get_edge_lonlats() + + np.testing.assert_allclose(lons, [-45., -90., -135., -135., -180., 135., + 135., 90., 45., 45., 0., -45.]) + np.testing.assert_allclose(lats, [80., 80., 80., 80., 80., 80., 80., + 80., 80., 80., 80., 80.]) + + def test_compute_optimal_bb(self): + """Test computing the bb area.""" + lats = np.array([[85.23900604248047, 62.256004333496094, 35.58000183105469], + [80.84000396728516, 60.74200439453125, 34.08500289916992], + [67.07600402832031, 54.147003173828125, 30.547000885009766]]).T + + lons = np.array([[-90.67900085449219, -21.565000534057617, -21.525001525878906], + [79.11000061035156, 7.284000396728516, -5.107000350952148], + [81.26400756835938, 29.672000885009766, 10.260000228881836]]).T + + area = geometry.SwathDefinition(lons, lats) + + res = area.compute_optimal_bb_area({'proj': 'omerc', 'ellps': 'WGS84'}) + + np.testing.assert_allclose(res.area_extent, (2286629.731529, + -2359693.817959, + 11729881.856072, + 2437001.523925)) + proj_dict = {'no_rot': True, 'lonc': 5.340645620216994, + 'ellps': 'WGS84', 'proj': 'omerc', + 'alpha': 19.022450179020247, 'lat_0': 60.7420043944989} + assert_np_dict_allclose(res.proj_dict, proj_dict) + self.assertEqual(res.shape, (3, 3)) + + class TestStackedAreaDefinition(unittest.TestCase): """Test the StackedAreaDefition.""" @@ -859,6 +959,86 @@ area_extent) +class TestDynamicAreaDefinition(unittest.TestCase): + """Test the DynamicAreaDefinition class.""" + + def test_freeze(self): + """Test freezing the area.""" + area = geometry.DynamicAreaDefinition('test_area', 'A test area', + {'proj': 'laea'}) + lons = [10, 10, 22, 22] + lats = [50, 66, 66, 50] + result = area.freeze((lons, lats), + resolution=3000, + proj_info={'lon0': 16, 'lat0': 58}) + + np.testing.assert_allclose(result.area_extent, (538546.7274949469, + 5380808.879250369, + 1724415.6519203288, + 6998895.701001488)) + self.assertEqual(result.proj_dict['lon0'], 16) + self.assertEqual(result.proj_dict['lat0'], 58) + self.assertEqual(result.x_size, 395) + self.assertEqual(result.y_size, 539) + + def test_freeze_with_bb(self): + """Test freezing the area with bounding box computation.""" + # area = geometry.DynamicAreaDefinition('test_area', 'A test area', + # {'proj': 'omerc'}, + # optimize_projection=False) + # lons = [[10, 12.1, 14.2, 16.3], + # [10, 12, 14, 16], + # [10, 11.9, 13.8, 15.7]] + # lats = [[66, 67, 68, 69.], + # [58, 59, 60, 61], + # [50, 51, 52, 53]] + # sdef = geometry.SwathDefinition(lons, lats) + # result = area.freeze(sdef, + # resolution=1000) + # self.assertTupleEqual(result.area_extent, (5578795.1654752363, + # -270848.61872542271, + # 7694893.3964453982, + # 126974.877141819)) + # self.assertEqual(result.x_size, 2116) + # self.assertEqual(result.y_size, 398) + + area = geometry.DynamicAreaDefinition('test_area', 'A test area', + {'proj': 'omerc'}, + optimize_projection=True) + lons = [[10, 12.1, 14.2, 16.3], + [10, 12, 14, 16], + [10, 11.9, 13.8, 15.7]] + lats = [[66, 67, 68, 69.], + [58, 59, 60, 61], + [50, 51, 52, 53]] + sdef = geometry.SwathDefinition(lons, lats) + result = area.freeze(sdef, + resolution=1000) + np.testing.assert_allclose(result.area_extent, (5050520.6077326955, + -336485.86803662963, + 8223167.9541879389, + 192612.12645302597)) + self.assertEqual(result.x_size, 3) + self.assertEqual(result.y_size, 4) + + def test_compute_domain(self): + """Test computing size and area extent.""" + area = geometry.DynamicAreaDefinition('test_area', 'A test area', + {'proj': 'laea'}) + corners = [1, 1, 9, 9] + self.assertRaises(ValueError, area.compute_domain, corners, 1, 1) + + area_extent, x_size, y_size = area.compute_domain(corners, size=(5, 5)) + self.assertTupleEqual(area_extent, (0, 0, 10, 10)) + self.assertEqual(x_size, 5) + self.assertEqual(y_size, 5) + + area_extent, x_size, y_size = area.compute_domain(corners, resolution=2) + self.assertTupleEqual(area_extent, (0, 0, 10, 10)) + self.assertEqual(x_size, 5) + self.assertEqual(y_size, 5) + + def suite(): """The test suite. """ @@ -866,6 +1046,8 @@ mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(Test)) mysuite.addTest(loader.loadTestsFromTestCase(TestStackedAreaDefinition)) + mysuite.addTest(loader.loadTestsFromTestCase(TestDynamicAreaDefinition)) + mysuite.addTest(loader.loadTestsFromTestCase(TestSwathDefinition)) return mysuite diff -Nru pyresample-1.5.0/pyresample/test/test_image.py pyresample-1.7.0/pyresample/test/test_image.py --- pyresample-1.5.0/pyresample/test/test_image.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/test_image.py 2017-10-13 09:24:57.000000000 +0000 @@ -18,7 +18,8 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 -Nru pyresample-1.5.0/pyresample/test/test_kd_tree.py pyresample-1.7.0/pyresample/test/test_kd_tree.py --- pyresample-1.5.0/pyresample/test/test_kd_tree.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/test_kd_tree.py 2017-10-13 09:24:57.000000000 +0000 @@ -4,8 +4,9 @@ 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 @@ 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 @@ 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') @@ -345,17 +346,19 @@ self.assertTrue(any(['Possible more' in str( x.message) for x in w]), 'Failed to create correct neighbour radius warning') cross_sum = res.sum() - cross_sum_stddev = stddev.sum() cross_sum_counts = counts.sum() expected = 1461.84313918 - expected_stddev = 0.446204424799 + expected_stddev = [0.446193170875, 0.443606880035, 0.438586349519] expected_counts = 4934802.0 self.assertTrue(res.shape == stddev.shape and stddev.shape == counts.shape and counts.shape == (800, 800, 3)) self.assertAlmostEqual(cross_sum, expected, msg='Swath multi channel resampling gauss failed on data') - self.assertAlmostEqual(cross_sum_stddev, expected_stddev, - msg='Swath multi channel resampling gauss failed on stddev') + for i, e_stddev in enumerate(expected_stddev): + cross_sum_stddev = stddev[:, :, i].sum() + print(cross_sum_stddev, e_stddev) + self.assertAlmostEqual(cross_sum_stddev, e_stddev, + msg='Swath multi channel resampling gauss failed on stddev (channel {})'.format(i)) self.assertAlmostEqual(cross_sum_counts, expected_counts, msg='Swath multi channel resampling gauss failed on counts') @@ -703,7 +706,7 @@ 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, @@ -822,3 +825,6 @@ mysuite.addTest(loader.loadTestsFromTestCase(Test)) return mysuite + +if __name__ == '__main__': + unittest.main() diff -Nru pyresample-1.5.0/pyresample/test/test_utils.py pyresample-1.7.0/pyresample/test/test_utils.py --- pyresample-1.5.0/pyresample/test/test_utils.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/test_utils.py 2017-10-13 09:24:57.000000000 +0000 @@ -3,7 +3,7 @@ import numpy as np -from pyresample import utils +from pyresample.test.utils import create_test_longitude, create_test_latitude def tmp(f): @@ -11,10 +11,10 @@ return f -class Test(unittest.TestCase): - +class TestLegacyAreaParser(unittest.TestCase): def test_area_parser_legacy(self): """Test legacy area parser.""" + from pyresample import utils ease_nh, ease_sh = utils.parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh', 'ease_sh') @@ -37,8 +37,32 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""" self.assertEquals(ease_sh.__str__(), sh_str) + def test_load_area(self): + from pyresample import utils + ease_nh = utils.load_area(os.path.join(os.path.dirname(__file__), + 'test_files', + 'areas.cfg'), 'ease_nh') + nh_str = """Area ID: ease_nh +Description: Arctic EASE grid +Projection ID: ease_nh +Projection: {'a': '6371228.0', 'lat_0': '90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'} +Number of columns: 425 +Number of rows: 425 +Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""" + self.assertEquals(nh_str, ease_nh.__str__()) + + def test_not_found_exception(self): + from pyresample import utils + self.assertRaises(utils.AreaNotFound, utils.parse_area_file, + os.path.join( + os.path.dirname(__file__), 'test_files', 'areas.cfg'), + 'no_area') + + +class TestYAMLAreaParser(unittest.TestCase): def test_area_parser_yaml(self): """Test YAML area parser.""" + from pyresample import utils ease_nh, ease_sh = utils.parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.yaml'), @@ -60,27 +84,102 @@ Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""" self.assertEquals(ease_sh.__str__(), sh_str) - def test_load_area(self): - ease_nh = utils.load_area(os.path.join(os.path.dirname(__file__), - 'test_files', - 'areas.cfg'), 'ease_nh') - nh_str = """Area ID: ease_nh -Description: Arctic EASE grid -Projection ID: ease_nh -Projection: {'a': '6371228.0', 'lat_0': '90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'} -Number of columns: 425 -Number of rows: 425 -Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""" - self.assertEquals(nh_str, ease_nh.__str__()) + def test_multiple_file_content(self): + from pyresample import utils + area_list = ["""ease_sh: + description: Antarctic EASE grid + projection: + a: 6371228.0 + units: m + lon_0: 0 + proj: laea + lat_0: -90 + shape: + height: 425 + width: 425 + area_extent: + lower_left_xy: [-5326849.0625, -5326849.0625] + upper_right_xy: [5326849.0625, 5326849.0625] + units: m +""", + """ease_sh2: + description: Antarctic EASE grid + projection: + a: 6371228.0 + units: m + lon_0: 0 + proj: laea + lat_0: -90 + shape: + height: 425 + width: 425 + area_extent: + lower_left_xy: [-5326849.0625, -5326849.0625] + upper_right_xy: [5326849.0625, 5326849.0625] + units: m +"""] + results = utils.parse_area_file(area_list) + self.assertEquals(len(results), 2) + self.assertIn(results[0].area_id, ('ease_sh', 'ease_sh2')) + self.assertIn(results[1].area_id, ('ease_sh', 'ease_sh2')) + + +class TestPreprocessing(unittest.TestCase): + def test_nearest_neighbor_area_area(self): + from pyresample import utils, geometry + proj_str = "+proj=lcc +datum=WGS84 +ellps=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs" + proj_dict = utils.proj4_str_to_dict(proj_str) + extents = [0, 0, 1000. * 5000, 1000. * 5000] + area_def = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS', + proj_dict, 400, 500, extents) + + extents2 = [-1000, -1000, 1000. * 4000, 1000. * 4000] + area_def2 = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS', + proj_dict, 600, 700, extents2) + rows, cols = utils.generate_nearest_neighbour_linesample_arrays(area_def, area_def2, 12000.) + + def test_nearest_neighbor_area_grid(self): + from pyresample import utils, geometry + lon_arr = create_test_longitude(-94.9, -90.0, (50, 100), dtype=np.float64) + lat_arr = create_test_latitude(25.1, 30.0, (50, 100), dtype=np.float64) + grid = geometry.GridDefinition(lons=lon_arr, lats=lat_arr) + + proj_str = "+proj=lcc +datum=WGS84 +ellps=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs" + proj_dict = utils.proj4_str_to_dict(proj_str) + extents = [0, 0, 1000. * 5000, 1000. * 5000] + area_def = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS', + proj_dict, 400, 500, extents) + rows, cols = utils.generate_nearest_neighbour_linesample_arrays(area_def, grid, 12000.) + + def test_nearest_neighbor_grid_area(self): + from pyresample import utils, geometry + proj_str = "+proj=lcc +datum=WGS84 +ellps=WGS84 +lat_0=25 +lat_1=25 +lon_0=-95 +units=m +no_defs" + proj_dict = utils.proj4_str_to_dict(proj_str) + extents = [0, 0, 1000. * 2500., 1000. * 2000.] + area_def = geometry.AreaDefinition('CONUS', 'CONUS', 'CONUS', + proj_dict, 40, 50, extents) + + lon_arr = create_test_longitude(-100.0, -60.0, (550, 500), dtype=np.float64) + lat_arr = create_test_latitude(20.0, 45.0, (550, 500), dtype=np.float64) + grid = geometry.GridDefinition(lons=lon_arr, lats=lat_arr) + rows, cols = utils.generate_nearest_neighbour_linesample_arrays(grid, area_def, 12000.) + + def test_nearest_neighbor_grid_grid(self): + from pyresample import utils, geometry + lon_arr = create_test_longitude(-95.0, -85.0, (40, 50), dtype=np.float64) + lat_arr = create_test_latitude(25.0, 35.0, (40, 50), dtype=np.float64) + grid_dst = geometry.GridDefinition(lons=lon_arr, lats=lat_arr) + + lon_arr = create_test_longitude(-100.0, -80.0, (400, 500), dtype=np.float64) + lat_arr = create_test_latitude(20.0, 40.0, (400, 500), dtype=np.float64) + grid = geometry.GridDefinition(lons=lon_arr, lats=lat_arr) + rows, cols = utils.generate_nearest_neighbour_linesample_arrays(grid, grid_dst, 12000.) - def test_not_found_exception(self): - self.assertRaises(utils.AreaNotFound, utils.parse_area_file, - os.path.join( - os.path.dirname(__file__), 'test_files', 'areas.cfg'), - 'no_area') +class TestMisc(unittest.TestCase): def test_wrap_longitudes(self): # test that we indeed wrap to [-180:+180[ + from pyresample import utils step = 60 lons = np.arange(-360, 360 + step, step) self.assertTrue( @@ -92,16 +191,53 @@ def test_unicode_proj4_string(self): """Test that unicode is accepted for area creation. """ + from pyresample import utils utils.get_area_def(u"eurol", u"eurol", u"bla", u'+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45', 1000, 1000, (-1000, -1000, 1000, 1000)) + def test_proj4_radius_parameters_provided(self): + from pyresample import utils + a, b = utils.proj4_radius_parameters( + '+proj=stere +a=6378273 +b=6356889.44891', + ) + np.testing.assert_almost_equal(a, 6378273) + np.testing.assert_almost_equal(b, 6356889.44891) + + def test_proj4_radius_parameters_ellps(self): + from pyresample import utils + a, b = utils.proj4_radius_parameters( + '+proj=stere +ellps=WGS84', + ) + np.testing.assert_almost_equal(a, 6378137.) + np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) + + def test_proj4_radius_parameters_default(self): + from pyresample import utils + a, b = utils.proj4_radius_parameters( + '+proj=lcc', + ) + # WGS84 + np.testing.assert_almost_equal(a, 6378137.) + np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) + + def 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. """ loader = unittest.TestLoader() mysuite = unittest.TestSuite() - mysuite.addTest(loader.loadTestsFromTestCase(Test)) + mysuite.addTest(loader.loadTestsFromTestCase(TestLegacyAreaParser)) + mysuite.addTest(loader.loadTestsFromTestCase(TestYAMLAreaParser)) + mysuite.addTest(loader.loadTestsFromTestCase(TestPreprocessing)) + mysuite.addTest(loader.loadTestsFromTestCase(TestMisc)) return mysuite diff -Nru pyresample-1.5.0/pyresample/test/utils.py pyresample-1.7.0/pyresample/test/utils.py --- pyresample-1.5.0/pyresample/test/utils.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/test/utils.py 2017-10-13 09:24:57.000000000 +0000 @@ -25,6 +25,7 @@ import types import warnings +import numpy as np _deprecations_as_exceptions = False _include_astropy_deprecations = False @@ -151,3 +152,24 @@ treat_deprecations_as_exceptions() +def create_test_longitude(start, stop, shape, twist_factor=0.0, dtype=np.float32): + if start > stop: + stop += 360.0 + + lon_row = np.linspace(start, stop, num=shape[1]).astype(dtype) + twist_array = np.arange(shape[0]).reshape((shape[0], 1)) * twist_factor + lon_array = np.repeat([lon_row], shape[0], axis=0) + lon_array += twist_array + + if stop > 360.0: + lon_array[lon_array > 360.0] -= 360 + return lon_array + + +def create_test_latitude(start, stop, shape, twist_factor=0.0, dtype=np.float32): + lat_col = np.linspace(start, stop, num=shape[0]).astype(dtype).reshape((shape[0], 1)) + twist_array = np.arange(shape[1]) * twist_factor + lat_array = np.repeat(lat_col, shape[1], axis=1) + lat_array += twist_array + return lat_array + diff -Nru pyresample-1.5.0/pyresample/utils.py pyresample-1.7.0/pyresample/utils.py --- pyresample-1.5.0/pyresample/utils.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/utils.py 2017-10-13 09:24:57.000000000 +0000 @@ -23,10 +23,12 @@ from __future__ import absolute_import +import os import numpy as np import six import yaml from configobj import ConfigObj +from collections import Mapping import pyresample as pr @@ -95,12 +97,36 @@ return _parse_legacy_area_file(area_file_name, *regions) -def _parse_yaml_area_file(area_file_name, *regions): - """Parse area information from a yaml area file.""" +def _read_yaml_area_file_content(area_file_name): + """Read one or more area files in to a single dict object.""" + if isinstance(area_file_name, (str, six.text_type)): + area_file_name = [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)): + 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 - with open(area_file_name) as fd_: - area_dict = yaml.load(fd_) +def _parse_yaml_area_file(area_file_name, *regions): + """Parse area information from a yaml area file. + + Args: + area_file_name: filename, file-like object, yaml string, or list of + these. + + The result of loading multiple area files is the combination of all + the files, using the first file as the "base", replacing things after + that. + """ + area_dict = _read_yaml_area_file_content(area_file_name) area_list = regions or area_dict.keys() res = [] @@ -113,21 +139,53 @@ area_name, area_file_name)) description = params['description'] projection = params['projection'] - xsize = params['shape']['width'] - ysize = params['shape']['height'] - area_extent = (params['area_extent']['lower_left_xy'] + - params['area_extent']['upper_right_xy']) - res.append(pr.geometry.AreaDefinition(area_name, description, - None, projection, - xsize, ysize, - area_extent)) + optimize_projection = params.get('optimize_projection', False) + try: + xsize = params['shape']['width'] + ysize = params['shape']['height'] + except KeyError: + xsize, ysize = None, None + try: + area_extent = (params['area_extent']['lower_left_xy'] + + params['area_extent']['upper_right_xy']) + except KeyError: + area_extent = None + area = pr.geometry.DynamicAreaDefinition(area_name, description, + projection, xsize, ysize, + area_extent, + optimize_projection) + try: + area = area.freeze() + except (TypeError, AttributeError): + pass + + res.append(area) return res +def _read_legacy_area_file_lines(area_file_name): + if isinstance(area_file_name, (str, six.text_type)): + area_file_name = [area_file_name] + + for area_file_obj in area_file_name: + if (isinstance(area_file_obj, (str, six.text_type)) and + not os.path.isfile(area_file_obj)): + # file content string + for line in area_file_obj.splitlines(): + yield line + continue + elif isinstance(area_file_obj, (str, six.text_type)): + # filename + with open(area_file_obj, 'r') as area_file_obj: + + for line in area_file_obj.readlines(): + yield line + + def _parse_legacy_area_file(area_file_name, *regions): """Parse area information from a legacy area file.""" - area_file = open(area_file_name, 'r') + area_file = _read_legacy_area_file_lines(area_file_name) area_list = list(regions) if len(area_list) == 0: select_all_areas = True @@ -138,7 +196,7 @@ # Extract area from file in_area = False - for line in area_file.readlines(): + for line in area_file: if not in_area: if 'REGION' in line: area_id = line.replace('REGION:', ''). \ @@ -156,8 +214,6 @@ else: area_content += line - area_file.close() - # Check if all specified areas were found if not select_all_areas: for i, area in enumerate(area_defs): @@ -227,19 +283,20 @@ """ 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 ----------- source_area_def : object - Source area definition as AreaDefinition object + Source area definition as geometry definition object target_area_def : object - Target area definition as AreaDefinition object + Target area definition as geometry definition object nprocs : int, optional Number of processor cores to be used @@ -247,11 +304,6 @@ ------- (row_indices, col_indices) : tuple of numpy arrays """ - if not (isinstance(source_area_def, pr.geometry.AreaDefinition) and - isinstance(target_area_def, pr.geometry.AreaDefinition)): - raise TypeError('source_area_def and target_area_def must be of type ' - 'geometry.AreaDefinition') - lons, lats = target_area_def.get_lonlats(nprocs) source_pixel_y, source_pixel_x = pr.grid.get_linesample(lons, lats, @@ -266,16 +318,18 @@ 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 ----------- source_area_def : object - Source area definition as AreaDefinition object + Source area definition as geometry definition object target_area_def : object - Target area definition as AreaDefinition object + Target area definition as geometry definition object radius_of_influence : float Cut off distance in meters nprocs : int, optional @@ -286,11 +340,6 @@ (row_indices, col_indices) : tuple of numpy arrays """ - if not (isinstance(source_area_def, pr.geometry.AreaDefinition) and - isinstance(target_area_def, pr.geometry.AreaDefinition)): - raise TypeError('source_area_def and target_area_def must be of type ' - 'geometry.AreaDefinition') - valid_input_index, valid_output_index, index_array, distance_array = \ pr.kd_tree.get_neighbour_info(source_area_def, target_area_def, @@ -359,13 +408,73 @@ def proj4_str_to_dict(proj4_str): """Convert PROJ.4 compatible string definition to dict - + Note: Key only parameters will be assigned a value of `True`. """ - pairs = (x.split('=', 1) for x in proj4_str.split(" ")) + 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. + + Arguments: + proj4_dict (str or dict): PROJ.4 parameters + + Returns: + a (float), b (float): equatorial and polar radius + """ + if isinstance(proj4_dict, str): + new_info = proj4_str_to_dict(proj4_dict) + else: + new_info = proj4_dict.copy() + + # load information from PROJ.4 about the ellipsis if possible + + 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: + 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): """Try to downcast array to uint16 """ @@ -393,3 +502,20 @@ """ lons_wrap = (lons + 180) % (360) - 180 return lons_wrap.astype(lons.dtype) + + +def recursive_dict_update(d, u): + """Recursive dictionary update using + + Copied from: + + http://stackoverflow.com/questions/3232943/update-value-of-a-nested-dictionary-of-varying-depth + + """ + for k, v in u.items(): + if isinstance(v, Mapping): + r = recursive_dict_update(d.get(k, {}), v) + d[k] = r + else: + d[k] = u[k] + return d diff -Nru pyresample-1.5.0/pyresample/version.py pyresample-1.7.0/pyresample/version.py --- pyresample-1.5.0/pyresample/version.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/pyresample/version.py 2017-10-13 09:24:57.000000000 +0000 @@ -15,4 +15,4 @@ # You should have received a copy of the GNU Lesser General Public License along # with this program. If not, see . -__version__ = '1.5.0' +__version__ = '1.7.0' diff -Nru pyresample-1.5.0/setup.cfg pyresample-1.7.0/setup.cfg --- pyresample-1.5.0/setup.cfg 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/setup.cfg 2017-10-13 09:24:57.000000000 +0000 @@ -1,5 +1,4 @@ [bdist_rpm] -requires=numpy pykdtree numexpr pyproj python-configobj +requires=numpy pykdtree python2-numexpr pyproj python-configobj release=1 doc_files = docs/Makefile docs/source/*.rst - diff -Nru pyresample-1.5.0/setup.py pyresample-1.7.0/setup.py --- pyresample-1.5.0/setup.py 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/setup.py 2017-10-13 09:24:57.000000000 +0000 @@ -27,11 +27,14 @@ version = imp.load_source('pyresample.version', 'pyresample/version.py') requirements = ['setuptools>=3.2', 'pyproj', 'numpy', 'configobj', - 'pykdtree>=1.1.1', 'pyyaml'] + 'pykdtree>=1.1.1', 'pyyaml', 'six'] extras_require = {'pykdtree': ['pykdtree>=1.1.1'], 'numexpr': ['numexpr'], 'quicklook': ['matplotlib', 'basemap', 'pillow']} +test_requires = [] +if sys.version_info < (3, 3): + test_requires.append('mock') if sys.version_info < (2, 6): # multiprocessing is not in the standard library requirements.append('multiprocessing') @@ -110,6 +113,7 @@ setup_requires=['numpy'], install_requires=requirements, extras_require=extras_require, + tests_require=test_requires, cmdclass={'build_ext': build_ext}, ext_modules=cythonize(extensions), test_suite='pyresample.test.suite', diff -Nru pyresample-1.5.0/.travis.yml pyresample-1.7.0/.travis.yml --- pyresample-1.5.0/.travis.yml 2017-05-02 08:24:36.000000000 +0000 +++ pyresample-1.7.0/.travis.yml 2017-10-13 09:24:57.000000000 +0000 @@ -1,28 +1,29 @@ language: python python: - '2.7' -- '3.3' - '3.4' - '3.5' +- '3.6' before_install: - sudo add-apt-repository ppa:ubuntugis/ppa -y - sudo apt-get update -qq - sudo apt-get install libfreetype6-dev -- sudo apt-get install libgeos-3.3.8 libgeos-c1 libgeos-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,!=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: