diff --git a/.travis.yml b/.travis.yml index 271de55..8ac7144 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,3 +1,4 @@ +dist: xenial language: python sudo: false notifications: @@ -9,8 +10,8 @@ addons: python: # We don't actually use the Travis Python, but this keeps it organized. - "2.7" - - "3.5" - "3.6" + - "3.7" install: # You may want to periodically update this, although the conda update # conda line below will keep everything up-to-date. We do this @@ -29,10 +30,11 @@ install: # Useful for debugging any issues with conda - conda info -a - - conda create -q -n test-environment -c conda-forge python=$TRAVIS_PYTHON_VERSION numpy scipy - - source activate test-environment - - pip install pytesmo pyresample - - pip install . + - conda create -q -n repurpose -c conda-forge python=$TRAVIS_PYTHON_VERSION + - conda env update -n repurpose -f environment.yml + - source activate repurpose + - python setup.py install + - pip list - which pip - which python diff --git a/CHANGES.rst b/CHANGES.rst index 7d41b65..03da615 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,6 +2,11 @@ Changelog ========= +Unreleased +========== + +- Add resample functions (from pytesmo) + Version 0.6 =========== diff --git a/README.rst b/README.rst index 7a61027..702fa69 100644 --- a/README.rst +++ b/README.rst @@ -60,6 +60,8 @@ It includes two main modules: spatial resampling. - ``ts2img`` for time series to image conversion, including support for temporal resampling. This module is very experimental at the moment. +- ``resample`` for spatial resampling of (regular or irregular) gridded data +to different resolutions. Alternatives ============ diff --git a/docs/index.rst b/docs/index.rst index 37458bc..45742ba 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,4 +1,6 @@ .. include:: ../README.rst +.. include:: ts2img.rst +.. include:: resample.rst Contents ======== diff --git a/docs/resample.rst b/docs/resample.rst new file mode 100644 index 0000000..f6fd263 --- /dev/null +++ b/docs/resample.rst @@ -0,0 +1,11 @@ +.. _resample: + +======== +resample +======== + +Spatial Resampling +================== + +The ``resample`` module contains functions to convert gridded data to different +spatial resolutions. \ No newline at end of file diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..a4bbd96 --- /dev/null +++ b/environment.yml @@ -0,0 +1,15 @@ +# List of conda and pip packages that are should be install when developing the packages +# To create the full conda environment: conda env create -f environment.yml +name: repurpose +channels: + - conda-forge + - defaults +dependencies: + - numpy + - scipy + - pip + - pip: + - pygeogrids + - pynetcf + - pyresample + - more_itertools diff --git a/repurpose/img2ts.py b/repurpose/img2ts.py index cd7036a..0b83faa 100644 --- a/repurpose/img2ts.py +++ b/repurpose/img2ts.py @@ -35,7 +35,7 @@ import pynetcf.time_series as nc import pygeogrids.grids as grids -import pytesmo.grid.resample as resamp +import repurpose.resample as resamp import numpy as np import os from datetime import datetime diff --git a/repurpose/resample.py b/repurpose/resample.py new file mode 100755 index 0000000..6311513 --- /dev/null +++ b/repurpose/resample.py @@ -0,0 +1,296 @@ +# Copyright (c) 2014,Vienna University of Technology, Department of Geodesy and Geoinformation +# All rights reserved. + +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# * Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# * Neither the name of the Vienna University of Technology, Department of Geodesy and Geoinformation nor the +# names of its contributors may be used to endorse or promote products +# derived from this software without specific prior written permission. + +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL VIENNA UNIVERSITY OF TECHNOLOGY, +# DEPARTMENT OF GEODESY AND GEOINFORMATION BE LIABLE FOR ANY +# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +''' +Created on Mar 25, 2014 + +@author: Christoph Paulik christoph.paulik@geo.tuwien.ac.at +''' + + +from pyresample import geometry, kd_tree +import numpy as np + + +def resample_to_grid_only_valid_return(input_data, src_lon, src_lat, target_lon, target_lat, + methods='nn', weight_funcs=None, + min_neighbours=1, search_rad=18000, neighbours=8, + fill_values=None): + """ + resamples data from dictionary of numpy arrays using pyresample + to given grid. + Searches for the neighbours and then resamples the data + to the grid given in togrid if at least + min_neighbours neighbours are found + + Parameters + ---------- + input_data : dict of numpy.arrays + src_lon : numpy.array + longitudes of the input data + src_lat : numpy.array + src_latitudes of the input data + target_lon : numpy.array + longitudes of the output data + target_src_lat : numpy.array + src_latitudes of the output data + methods : string or dict, optional + method of spatial averaging. this is given to pyresample + and can be + 'nn' : nearest neighbour + 'custom' : custom weight function has to be supplied in weight_funcs + see pyresample documentation for more details + can also be a dictionary with a method for each array in input data dict + weight_funcs : function or dict of functions, optional + if method is 'custom' a function like func(distance) has to be given + can also be a dictionary with a function for each array in input data dict + min_neighbours: int, optional + if given then only points with at least this number of neighbours will be + resampled + Default : 1 + search_rad : float, optional + search radius in meters of neighbour search + Default : 18000 + neighbours : int, optional + maximum number of neighbours to look for for each input grid point + Default : 8 + fill_values : number or dict, optional + if given the output array will be filled with this value if no valid + resampled value could be computed, if not a masked array will be returned + can also be a dict with a fill value for each variable + Returns + ------- + data : dict of numpy.arrays + resampled data on part of the target grid over which data was found + mask: numpy.ndarray + boolean mask into target grid that specifies where data was resampled + + Raises + ------ + ValueError : + if empty dataset is resampled + """ + output_data = {} + + if target_lon.ndim == 2: + target_lat = target_lat.ravel() + target_lon = target_lon.ravel() + + input_swath = geometry.SwathDefinition(src_lon, src_lat) + output_swath = geometry.SwathDefinition(target_lon, target_lat) + + (valid_input_index, + valid_output_index, + index_array, + distance_array) = kd_tree.get_neighbour_info(input_swath, + output_swath, + search_rad, + neighbours=neighbours) + + # throw away points with less than min_neighbours neighbours + # find points with valid neighbours + # get number of found neighbours for each grid point/row + if neighbours > 1: + nr_neighbours = np.isfinite(distance_array).sum(1) + neigh_condition = nr_neighbours >= min_neighbours + mask = np.invert(neigh_condition) + enough_neighbours = np.nonzero(neigh_condition)[0] + if neighbours == 1: + nr_neighbours = np.isfinite(distance_array) + neigh_condition = nr_neighbours >= min_neighbours + mask = np.invert(neigh_condition) + enough_neighbours = np.nonzero(neigh_condition)[0] + distance_array = np.reshape( + distance_array, (distance_array.shape[0], 1)) + index_array = np.reshape(index_array, (index_array.shape[0], 1)) + + if enough_neighbours.size == 0: + raise ValueError( + "No points with at least %d neighbours found" % min_neighbours) + + # remove neighbourhood info of input grid points that have no neighbours to not have to + # resample to whole output grid for small input grid file + distance_array = distance_array[enough_neighbours, :] + index_array = index_array[enough_neighbours, :] + valid_output_index = valid_output_index[enough_neighbours] + + for param in input_data: + + data = input_data[param] + + if type(methods) == dict: + method = methods[param] + else: + method = methods + + if method is not 'nn': + if type(weight_funcs) == dict: + weight_func = weight_funcs[param] + else: + weight_func = weight_funcs + else: + weight_func = None + + neigh_slice = slice(None, None, None) + # check if method is nn, if so only use first row of index_array and + # distance_array + if method == 'nn': + neigh_slice = (slice(None, None, None), 0) + + if type(fill_values) == dict: + fill_value = fill_values[param] + else: + fill_value = fill_values + + output_array = kd_tree.get_sample_from_neighbour_info( + method, + enough_neighbours.shape, + data, + valid_input_index, + valid_output_index, + index_array[neigh_slice], + distance_array[neigh_slice], + weight_funcs=weight_func, + fill_value=fill_value) + + output_data[param] = output_array + + return output_data, mask + + +def resample_to_grid(input_data, src_lon, src_lat, target_lon, target_lat, + methods='nn', weight_funcs=None, + min_neighbours=1, search_rad=18000, neighbours=8, + fill_values=None): + """ + resamples data from dictionary of numpy arrays using pyresample + to given grid. + Searches for the neighbours and then resamples the data + to the grid given in togrid if at least + min_neighbours neighbours are found + + Parameters + ---------- + input_data : dict of numpy.arrays + src_lon : numpy.array + longitudes of the input data + src_lat : numpy.array + src_latitudes of the input data + target_lon : numpy.array + longitudes of the output data + target_src_lat : numpy.array + src_latitudes of the output data + methods : string or dict, optional + method of spatial averaging. this is given to pyresample + and can be + 'nn' : nearest neighbour + 'custom' : custom weight function has to be supplied in weight_funcs + see pyresample documentation for more details + can also be a dictionary with a method for each array in input data dict + weight_funcs : function or dict of functions, optional + if method is 'custom' a function like func(distance) has to be given + can also be a dictionary with a function for each array in input data dict + min_neighbours: int, optional + if given then only points with at least this number of neighbours will be + resampled + Default : 1 + search_rad : float, optional + search radius in meters of neighbour search + Default : 18000 + neighbours : int, optional + maximum number of neighbours to look for for each input grid point + Default : 8 + fill_values : number or dict, optional + if given the output array will be filled with this value if no valid + resampled value could be computed, if not a masked array will be returned + can also be a dict with a fill value for each variable + Returns + ------- + data : dict of numpy.arrays + resampled data on given grid + Raises + ------ + ValueError : + if empty dataset is resampled + """ + + output_data = {} + output_shape = target_lat.shape + if target_lon.ndim == 2: + target_lat = target_lat.ravel() + target_lon = target_lon.ravel() + + resampled_data, mask = resample_to_grid_only_valid_return(input_data, + src_lon, src_lat, + target_lon, target_lat, + methods=methods, + weight_funcs=weight_funcs, + min_neighbours=min_neighbours, + search_rad=search_rad, + neighbours=neighbours) + for param in input_data: + data = resampled_data[param] + orig_data = input_data[param] + + if type(fill_values) == dict: + fill_value = fill_values[param] + else: + fill_value = fill_values + + # construct arrays in output grid form + if fill_value is not None: + output_array = np.zeros( + target_lat.shape, dtype=orig_data.dtype) + fill_value + else: + output_array = np.zeros(target_lat.shape, dtype=orig_data.dtype) + output_array = np.ma.array(output_array, mask=mask) + output_array[~mask] = data + + output_data[param] = output_array.reshape(output_shape) + + return output_data + + +def hamming_window(radius, distances): + """ + Hamming window filter. + + Parameters + ---------- + radius : float32 + Radius of the window. + distances : numpy.ndarray + Array with distances. + + Returns + ------- + weights : numpy.ndarray + Distance weights. + """ + alpha = 0.54 + weights = alpha + (1 - alpha) * np.cos(np.pi / radius * distances) + + return weights diff --git a/requirements.txt b/requirements.txt index d8a249d..b1a509f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ # Add your requirements here like: numpy pygeogrids>=0.1.3 -pytesmo>=0.3.0 pynetcf +pyresample +more_itertools diff --git a/test-requirements.txt b/test-requirements.txt index 468f195..c55812c 100644 --- a/test-requirements.txt +++ b/test-requirements.txt @@ -2,4 +2,5 @@ # They will be installed automatically when running `python setup.py test`. # ATTENTION: Don't remove pytest-cov and pytest as they are needed. pytest-cov -pytest +coverage==4.5.2 +pytest==3.8.2 diff --git a/test_requirements.txt b/test_requirements.txt deleted file mode 100644 index 4023b81..0000000 --- a/test_requirements.txt +++ /dev/null @@ -1,5 +0,0 @@ -# requirements for testing the package -# I think it makes little sense to mock the writer classes in the tests -# so we will use real classes -pynetcf -pygeobase diff --git a/tests/test_resample.py b/tests/test_resample.py new file mode 100755 index 0000000..a9b9011 --- /dev/null +++ b/tests/test_resample.py @@ -0,0 +1,137 @@ +# Copyright (c) 2014,Vienna University of Technology, Department of Geodesy and Geoinformation +# All rights reserved. + +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# * Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# * Neither the name of the Vienna University of Technology, Department of +# Geodesy and Geoinformation nor the +# names of its contributors may be used to endorse or promote products +# derived from this software without specific prior written permission. + +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL VIENNA UNIVERSITY OF TECHNOLOGY, +# DEPARTMENT OF GEODESY AND GEOINFORMATION BE LIABLE FOR ANY +# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +''' +Created on Jun 11, 2014 + +@author: Christoph Paulik christoph.paulik@geo.tuwien.ac.at +''' + + +import repurpose.resample as resample +import numpy as np +import functools +import numpy.testing as nptest +import unittest + + +class Test_resample(unittest.TestCase): + + def setUp(self): + nlats = 90*2 + nlons = 180*2 + self.testdata = dict(valrange=np.arange(0.,nlons*nlats).reshape(nlats,nlons), + ones=np.ones((nlats,nlons))) + self.lons, self.lats = np.meshgrid(np.arange(-180., 180., 1.), + np.arange(-90., 90., 1.)) + + def tearDown(self): + self.testdata = None + + def test_resample_to_zero_dot_one_deg(self): + # lets resample to a 0.5 degree grid + # define the grid points in latitude and longitude + lats_dim = np.arange(0, 0.4, 0.1) + lons_dim = np.arange(0, 0.4, 0.1) + # make 2d grid out the 1D grid spacing + lons_grid, lats_grid = np.meshgrid(lons_dim, lats_dim) + + resampled_data = resample.resample_to_grid(self.testdata, self.lons, self.lats, + lons_grid, lats_grid, search_rad=60000, + neighbours=1,fill_values=np.nan) + + for key in self.testdata: + assert resampled_data[key].shape == lons_grid.shape + + assert np.all(np.all(resampled_data['ones'], 1)) + assert np.all(resampled_data['valrange'] == self.testdata['valrange'][90, 180]) + + +def test_resample_dtypes(): + """ + Test if dtypes stay the same when resampling. + """ + + data = {'testint8': np.array([5, 5], dtype=np.int8), + 'testfloat16': np.array([5, 5], dtype=np.float16)} + + fill_values = {'testint8': 0, + 'testfloat16': 999.} + lons = np.array([0, 0.1]) + lats = np.array([0, 0.1]) + # lets resample to a 0.1 degree grid + # define the grid points in latitude and longitude + lats_dim = np.arange(-1, 1, 0.1) + lons_dim = np.arange(-1, 1, 0.1) + # make 2d grid out the 1D grid spacing + lons_grid, lats_grid = np.meshgrid(lons_dim, lats_dim) + + resampled_data = resample.resample_to_grid(data, lons, lats, + lons_grid, lats_grid, + fill_values=fill_values) + + for key in data: + assert resampled_data[key].shape == lons_grid.shape + assert resampled_data[key].dtype == data[key].dtype + + +def test_resample_hamming(): + """ + Test if hamming window is applied correctly + """ + # let's do 5 points with the highest value in the middle + # -1- + # 151 + # -1- + + data = {'testfloat16': np.array([1, 1, 5, 1, 1], dtype=np.float16)} + + fill_values = {'testfloat16': 999.} + lons = np.array([0, -0.1, 0, 0.1, 0]) + lats = np.array([0.1, 0, 0, 0, -0.1]) + # lets resample to a 0.1 degree grid + # define the grid points in latitude and longitude + lats_dim = np.arange(-0.1, 0.11, 0.1) + lons_dim = np.arange(-0.1, 0.11, 0.1) + # make 2d grid out the 1D grid spacing + lons_grid, lats_grid = np.meshgrid(lons_dim, lats_dim) + # make partial function of the hamming window the radius of the hamming + # window is in meters not in degrees + hamm = functools.partial(resample.hamming_window, 15000) + + resampled_data = resample.resample_to_grid(data, lons, lats, + lons_grid, lats_grid, + fill_values=fill_values, + methods='custom', weight_funcs=hamm) + + resampled_should = np.array([[1.640625, 1.64160156, 1.640625], + [1.64160156, 3.11132812, 1.64160156], + [1.640625, 1.64160156, 1.640625]]) + for key in data: + assert resampled_data[key].shape == lons_grid.shape + assert resampled_data[key].dtype == data[key].dtype + nptest.assert_almost_equal(resampled_data[key], resampled_should)