[Git][debian-gis-team/python-geotiepoints][master] 6 commits: New upstream version 1.3.0
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Sat Sep 18 10:22:59 BST 2021
Antonio Valentino pushed to branch master at Debian GIS Project / python-geotiepoints
Commits:
c8251758 by Antonio Valentino at 2021-09-18T07:40:00+00:00
New upstream version 1.3.0
- - - - -
7b84c9e2 by Antonio Valentino at 2021-09-18T07:40:40+00:00
Update upstream source from tag 'upstream/1.3.0'
Update to upstream version '1.3.0'
with Debian dir 395e840f289e0576d82103ac2d7459403adb776a
- - - - -
c079eff4 by Antonio Valentino at 2021-09-18T07:48:04+00:00
New upstream release
- - - - -
a3bd4a1c by Antonio Valentino at 2021-09-18T07:50:19+00:00
Update d/copyright
- - - - -
d778c0a8 by Antonio Valentino at 2021-09-18T08:58:43+00:00
Update and rename patch
- - - - -
e5d72eb1 by Antonio Valentino at 2021-09-18T08:59:13+00:00
Set distribution to unstable
- - - - -
17 changed files:
- .github/PULL_REQUEST_TEMPLATE.md
- CHANGELOG.md
- README.md
- RELEASING.md
- debian/changelog
- debian/copyright
- debian/patches/0001-Fix-tests.patch → debian/patches/0001-Skip-tests-using-external-data.patch
- debian/patches/series
- geotiepoints/geointerpolator.py
- geotiepoints/interpolator.py
- geotiepoints/modisinterpolator.py
- + geotiepoints/simple_modis_interpolator.py
- + geotiepoints/tests/test_simple_modis_interpolator.py
- + geotiepoints/tests/utils.py
- geotiepoints/version.py
- setup.cfg
- setup.py
Changes:
=====================================
.github/PULL_REQUEST_TEMPLATE.md
=====================================
@@ -3,5 +3,5 @@
- [ ] Closes #xxxx <!-- remove if there is no corresponding issue, which should only be the case for minor changes -->
- [ ] Tests added <!-- for all bug fixes or enhancements -->
- [ ] Tests passed <!-- for all non-documentation changes) -->
- - [ ] Passes ``git diff origin/master **/*py | flake8 --diff`` <!-- remove if you did not edit any Python files -->
+ - [ ] Passes ``git diff origin/main **/*py | flake8 --diff`` <!-- remove if you did not edit any Python files -->
- [ ] Fully documented <!-- remove if this change should not be visible to users, e.g., if it is an internal clean-up, or if this is part of a larger project that will be documented later -->
=====================================
CHANGELOG.md
=====================================
@@ -1,3 +1,14 @@
+## Version 1.3.0 (2021/09/12)
+
+### Pull Requests Merged
+
+#### Features added
+
+* [PR 31](https://github.com/pytroll/python-geotiepoints/pull/31) - Add simple lon/lat based MODIS interpolation
+
+In this release 1 pull request was closed.
+
+
## Version 1.2.1 (2021/03/08)
### Issues Closed
=====================================
README.md
=====================================
@@ -1,9 +1,9 @@
python-geotiepoints
===================
-[](https://github.com/pytroll/python-geotiepoints/actions?query=workflow%3A%22CI%22)
-[](https://coveralls.io/github/pytroll/python-geotiepoints?branch=master)
-[](https://landscape.io/github/pytroll/python-geotiepoints/master)
+[](https://travis-ci.org/pytroll/python-geotiepoints)
+[](https://coveralls.io/github/pytroll/python-geotiepoints?branch=main)
+[](https://landscape.io/github/pytroll/python-geotiepoints/main)
Python-geotiepoints is a python module that interpolates (and extrapolates if
=====================================
RELEASING.md
=====================================
@@ -2,13 +2,13 @@
prerequisites: `pip install loghub setuptools twine`
-1. checkout master
+1. checkout main branch
2. pull from repo
3. run the unittests
4. run `loghub` and update the `CHANGELOG.md` file:
```
-loghub pytroll/python-geotiepoints --token $LOGHUB_GITHUB_TOKEN -st v0.8.0 -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes -plg backwards-incompatibility "Backward incompatible changes" -plg refactor "Refactoring"
+loghub pytroll/python-geotiepoints --token $LOGHUB_GITHUB_TOKEN -st v0.8.0 -plg bug "Bugs fixed" -plg enhancement "Features added" -plg documentation "Documentation changes" -plg backwards-incompatibility "Backward incompatible changes" -plg refactor "Refactoring"
```
This uses a `LOGHUB_GITHUB_TOKEN` environment variable. This must be created
=====================================
debian/changelog
=====================================
@@ -1,9 +1,16 @@
-python-geotiepoints (1.2.1-2) UNRELEASED; urgency=medium
+python-geotiepoints (1.3.0-1) unstable; urgency=medium
- * Team upload.
+ [ Bas Couwenberg ]
* Bump Standards-Version to 4.6.0, no changes.
- -- Bas Couwenberg <sebastic at debian.org> Wed, 08 Sep 2021 17:42:45 +0200
+ [ Antonio Valentino ]
+ * New upstream release.
+ * Update debian/copyright file.
+ * debian/patches:
+ - update 0001-Fix-tests.patch and rename it into
+ 0001-Skip-tests-using-external-data.patch.
+
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Sat, 18 Sep 2021 08:58:56 +0000
python-geotiepoints (1.2.1-1) unstable; urgency=medium
=====================================
debian/copyright
=====================================
@@ -17,7 +17,7 @@ Copyright: 2018, Brian Warner
License: CC0-1.0
Files: debian/*
-Copyright: 2018-2020 Antonio Valentino <antonio.valentino at tiscali.it>
+Copyright: 2018-2021 Antonio Valentino <antonio.valentino at tiscali.it>
License: GPL-3+
License: GPL-3+
=====================================
debian/patches/0001-Fix-tests.patch → debian/patches/0001-Skip-tests-using-external-data.patch
=====================================
@@ -1,13 +1,14 @@
From: Antonio Valentino <antonio.valentino at tiscali.it>
Date: Sun, 23 Dec 2018 10:36:27 +0000
-Subject: Fix tests
+Subject: Skip tests using external data
Forwarded: not-needed
---
- geotiepoints/tests/test_modis.py | 4 ++++
- geotiepoints/tests/test_modisinterpolator.py | 1 +
- setup.py | 2 +-
- 3 files changed, 6 insertions(+), 1 deletion(-)
+ geotiepoints/tests/test_modis.py | 4 ++++
+ geotiepoints/tests/test_modisinterpolator.py | 1 +
+ geotiepoints/tests/test_simple_modis_interpolator.py | 2 ++
+ setup.py | 2 +-
+ 4 files changed, 8 insertions(+), 1 deletion(-)
diff --git a/geotiepoints/tests/test_modis.py b/geotiepoints/tests/test_modis.py
index ef9c562..a2a1eac 100644
@@ -43,8 +44,28 @@ index f688dc7..57503a9 100644
class TestModisInterpolator(unittest.TestCase):
def test_modis(self):
h5f = h5py.File(FILENAME_DATA, 'r')
+diff --git a/geotiepoints/tests/test_simple_modis_interpolator.py b/geotiepoints/tests/test_simple_modis_interpolator.py
+index bbe340f..de382b3 100644
+--- a/geotiepoints/tests/test_simple_modis_interpolator.py
++++ b/geotiepoints/tests/test_simple_modis_interpolator.py
+@@ -75,6 +75,7 @@ def _load_250m_lonlat_expected_as_xarray_dask():
+ return lon250, lat250
+
+
++ at pytest.mark.skipif(not os.path.exists(FILENAME_DATA), reason="data file not available")
+ @pytest.mark.parametrize(
+ ("input_func", "exp_func", "interp_func"),
+ [
+@@ -105,6 +106,7 @@ def test_basic_interp(input_func, exp_func, interp_func):
+ assert not np.any(np.isnan(lats))
+
+
++ at pytest.mark.skipif(not os.path.exists(FILENAME_DATA), reason="data file not available")
+ def test_nonstandard_scan_size():
+ lon1, lat1 = _load_1km_lonlat_as_xarray_dask()
+ # remove 1 row from the end
diff --git a/setup.py b/setup.py
-index f206095..5af8e0c 100644
+index 85a5dc9..bce938b 100644
--- a/setup.py
+++ b/setup.py
@@ -64,7 +64,7 @@ if __name__ == "__main__":
=====================================
debian/patches/series
=====================================
@@ -1 +1 @@
-0001-Fix-tests.patch
+0001-Skip-tests-using-external-data.patch
=====================================
geotiepoints/geointerpolator.py
=====================================
@@ -1,42 +1,37 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
-
-# Copyright (c) 2013 Martin Raspaud
-
-# Author(s):
-
-# Martin Raspaud <martin.raspaud at smhi.se>
-
+#
+# Copyright (c) 2013-2021 Python-geotiepoints developers
+#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
-
+#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
-
+#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-"""Geographical interpolation (lon/lats).
-"""
+"""Geographical interpolation (lon/lats)."""
import numpy as np
-from numpy import arccos, sign, rad2deg, sqrt, arcsin
from geotiepoints.interpolator import Interpolator
EARTH_RADIUS = 6370997.0
+
class GeoInterpolator(Interpolator):
- """
- Handles interpolation of geolocation from a grid of tie points. It is
+ """Handles interpolation of geolocation from a grid of tie points.
+
+ It is
preferable to have tie-points out till the edges if the tiepoint grid, but
a method is provided to extrapolate linearly the tiepoints to the borders
of the grid. The extrapolation is done automatically if it seems necessary.
- Uses numpy, scipy, and optionally pyresample
+ Uses numpy, scipy, and optionally pyresample.
The constructor takes in the tiepointed data as *data*, the
*tiepoint_grid* and the desired *final_grid*. As optional arguments, one
@@ -45,62 +40,53 @@ class GeoInterpolator(Interpolator):
along the y axis (this affects how the extrapolator behaves). If
*chunksize* is set, don't forget to adjust the interpolation orders
accordingly: the interpolation is indeed done globaly (not chunkwise).
+
"""
- def __init__(self, lon_lat_data, *args, **kwargs):
- Interpolator.__init__(self, None, *args, **kwargs)
- self.lon_tiepoint = None
- self.lat_tiepoint = None
+ def __init__(self, lon_lat_data, *args, **kwargs):
try:
# Maybe it's a pyresample object ?
- self.set_tiepoints(lon_lat_data.lons, lon_lat_data.lats)
+ self.lon_tiepoint = lon_lat_data.lons
+ self.lat_tiepoint = lon_lat_data.lats
xyz = lon_lat_data.get_cartesian_coords()
- self.tie_data = [xyz[:, :, 0], xyz[:, :, 1], xyz[:, :, 2]]
-
+ tie_data = [xyz[:, :, 0], xyz[:, :, 1], xyz[:, :, 2]]
except AttributeError:
- self.set_tiepoints(lon_lat_data[0], lon_lat_data[1])
- lons_rad = np.radians(self.lon_tiepoint)
- lats_rad = np.radians(self.lat_tiepoint)
- x__ = EARTH_RADIUS * np.cos(lats_rad) * np.cos(lons_rad)
- y__ = EARTH_RADIUS * np.cos(lats_rad) * np.sin(lons_rad)
- z__ = EARTH_RADIUS * np.sin(lats_rad)
- self.tie_data = [x__, y__, z__]
-
+ self.lon_tiepoint = lon_lat_data[0]
+ self.lat_tiepoint = lon_lat_data[1]
+ x__, y__, z__ = lonlat2xyz(self.lon_tiepoint, self.lat_tiepoint)
+ tie_data = [x__, y__, z__]
- self.new_data = []
- for num in range(len(self.tie_data)):
- self.new_data.append([])
-
-
- def set_tiepoints(self, lon, lat):
- """Defines the lon,lat tie points.
- """
- self.lon_tiepoint = lon
- self.lat_tiepoint = lat
+ super().__init__(tie_data, *args, **kwargs)
def interpolate(self):
- newx, newy, newz = Interpolator.interpolate(self)
- lon = get_lons_from_cartesian(newx, newy)
- lat = get_lats_from_cartesian(newx, newy, newz)
+ """Run the interpolation."""
+ newx, newy, newz = super().interpolate()
+ lon, lat = xyz2lonlat(newx, newy, newz)
return lon, lat
-def get_lons_from_cartesian(x__, y__):
- """Get longitudes from cartesian coordinates.
- """
- return rad2deg(arccos(x__ / sqrt(x__ ** 2 + y__ ** 2))) * sign(y__)
-
-def get_lats_from_cartesian(x__, y__, z__, thr=0.8):
- """Get latitudes from cartesian coordinates.
- """
- # if we are at low latitudes - small z, then get the
- # latitudes only from z. If we are at high latitudes (close to the poles)
- # then derive the latitude using x and y:
-
- lats = np.where(np.logical_and(np.less(z__, thr * EARTH_RADIUS),
- np.greater(z__, -1. * thr * EARTH_RADIUS)),
- 90 - rad2deg(arccos(z__/EARTH_RADIUS)),
- sign(z__) *
- (90 - rad2deg(arcsin(sqrt(x__ ** 2 + y__ ** 2)
- / EARTH_RADIUS))))
- return lats
+def lonlat2xyz(lons, lats, radius=EARTH_RADIUS):
+ """Convert lons and lats to cartesian coordinates."""
+ lons_rad = np.deg2rad(lons)
+ lats_rad = np.deg2rad(lats)
+ x_coords = radius * np.cos(lats_rad) * np.cos(lons_rad)
+ y_coords = radius * np.cos(lats_rad) * np.sin(lons_rad)
+ z_coords = radius * np.sin(lats_rad)
+ return x_coords, y_coords, z_coords
+
+
+def xyz2lonlat(x__, y__, z__, radius=EARTH_RADIUS, thr=0.8, low_lat_z=True):
+ """Get longitudes from cartesian coordinates."""
+ lons = np.rad2deg(np.arccos(x__ / np.sqrt(x__ ** 2 + y__ ** 2))) * np.sign(y__)
+ lats = np.sign(z__) * (90 - np.rad2deg(np.arcsin(np.sqrt(x__ ** 2 + y__ ** 2) / radius)))
+ if low_lat_z:
+ # if we are at low latitudes - small z, then get the
+ # latitudes only from z. If we are at high latitudes (close to the poles)
+ # then derive the latitude using x and y:
+ lat_mask_cond = np.logical_and(
+ np.less(z__, thr * radius),
+ np.greater(z__, -1. * thr * radius))
+ lat_z_only = 90 - np.rad2deg(np.arccos(z__ / radius))
+ lats = np.where(lat_mask_cond, lat_z_only, lats)
+
+ return lons, lats
=====================================
geotiepoints/interpolator.py
=====================================
@@ -104,10 +104,7 @@ class Interpolator(object):
else:
self.tie_data = list(data)
- self.new_data = []
- for num in range(len(self.tie_data)):
- self.new_data.append([])
-
+ self.new_data = [[] for _ in self.tie_data]
self.kx_, self.ky_ = kx_, ky_
def fill_borders(self, *args):
=====================================
geotiepoints/modisinterpolator.py
=====================================
@@ -31,6 +31,8 @@ import dask.array as da
import numpy as np
import warnings
+from .geointerpolator import lonlat2xyz, xyz2lonlat
+
R = 6371.
# Aqua scan width and altitude in km
scan_width = 10.00017
@@ -199,19 +201,8 @@ class ModisInterpolator(object):
a_scan = (s_s + s_s * (1 - s_s) * c_exp_full + s_t * (1 - s_t) * c_ali_full)
res = []
-
- sublat = lat1[::16, ::16]
- sublon = lon1[::16, ::16]
- to_cart = abs(sublat).max() > 60 or (sublon.max() - sublon.min()) > 180
-
- if to_cart:
- datasets = lonlat2xyz(lon1, lat1)
- else:
- datasets = [lon1, lat1]
-
+ datasets = lonlat2xyz(lon1, lat1)
for data in datasets:
- data_attrs = data.attrs
- dims = data.dims
data = data.data
data = data.reshape((-1, cscan_len, cscan_full_width))
data_a, data_b, data_c, data_d = get_corners(data)
@@ -224,12 +215,9 @@ class ModisInterpolator(object):
data_2 = (1 - a_scan) * data_d + a_scan * data_c
data = (1 - a_track) * data_1 + a_track * data_2
- res.append(xr.DataArray(data, attrs=data_attrs, dims=dims))
-
- if to_cart:
- return xyz2lonlat(*res)
- else:
- return res
+ res.append(data)
+ lon, lat = xyz2lonlat(*res)
+ return xr.DataArray(lon, dims=lon1.dims), xr.DataArray(lat, dims=lat1.dims)
def modis_1km_to_250m(lon1, lat1, satz1):
@@ -264,21 +252,3 @@ def modis_5km_to_250m(lon1, lat1, satz1):
"may result in poor quality")
interp = ModisInterpolator(5000, 250, lon1.shape[1])
return interp.interpolate(lon1, lat1, satz1)
-
-
-def lonlat2xyz(lons, lats):
- """Convert lons and lats to cartesian coordinates."""
- 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 x_coords, y_coords, z_coords
-
-def xyz2lonlat(x__, y__, z__):
- """Get longitudes from cartesian coordinates.
- """
- R = 6370997.0
- lons = da.rad2deg(da.arccos(x__ / da.sqrt(x__ ** 2 + y__ ** 2))) * da.sign(y__)
- lats = da.sign(z__) * (90 - da.rad2deg(da.arcsin(da.sqrt(x__ ** 2 + y__ ** 2) / R)))
-
- return lons, lats
=====================================
geotiepoints/simple_modis_interpolator.py
=====================================
@@ -0,0 +1,222 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Copyright (c) 2021 Python-geotiepoints developers
+#
+# This file is part of python-geotiepoints.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""Interpolate MODIS 1km navigation arrays to 250m and 500m resolutions.
+
+The code used here is a rewrite of the IDL function ``MODIS_GEO_INTERP_250``
+used by Liam Gumley. It has been modified to convert coordinates to cartesian
+(X, Y, Z) coordinates first to avoid problems with the anti-meridian and poles.
+This code was originally part of the CSPP Polar2Grid project, but has been
+moved here for integration with Satpy and newer versions of Polar2Grid.
+
+This algorithm differs from the one in ``modisinterpolator`` as it only
+requires the original longitude and latitude arrays. This is useful in the
+case of reading the 250m or 500m MODIS L1b files or any MODIS L2 files without
+including the MOD03 geolocation file as there is no SensorZenith angle dataset
+in these files.
+
+"""
+
+from functools import wraps
+
+import numpy as np
+from scipy.ndimage.interpolation import map_coordinates
+from .geointerpolator import lonlat2xyz, xyz2lonlat
+
+try:
+ import dask.array as da
+except ImportError:
+ # if dask can't be imported then we aren't going to be given dask arrays
+ da = None
+
+try:
+ import xarray as xr
+except ImportError:
+ xr = None
+
+# MODIS has 10 rows of data in the array for every scan line
+ROWS_PER_SCAN = 10
+
+
+def scanline_mapblocks(func):
+ """Convert dask array inputs to appropriate map_blocks calls.
+
+ This function, applied as a decorator, will call the wrapped function
+ using dask's ``map_blocks``. It will rechunk inputs when necessary to make
+ sure that the input chunks are entire scanlines to avoid incorrect
+ interpolation.
+
+ """
+ @wraps(func)
+ def _wrapper(lon_data, lat_data, res_factor=4):
+ if lon_data.ndim != 2 or lat_data.ndim != 2:
+ raise ValueError("Expected 2D lon/lat arrays.")
+ if hasattr(lon_data, "compute"):
+ # assume it is dask or xarray with dask, ensure proper chunk size
+ # if DataArray get just the dask array
+ lon_dask = lon_data.data if hasattr(lon_data, "dims") else lon_data
+ lat_dask = lat_data.data if hasattr(lat_data, "dims") else lat_data
+ lon_dask, lat_dask = _rechunk_lonlat_if_needed(lon_dask, lat_dask)
+ new_lons, new_lats = _call_map_blocks_interp(func, lon_dask, lat_dask, res_factor)
+ if hasattr(lon_data, "dims"):
+ # recreate DataArrays
+ new_lons = xr.DataArray(new_lons, dims=lon_data.dims)
+ new_lats = xr.DataArray(new_lats, dims=lon_data.dims)
+ return new_lons, new_lats
+
+ return func(lon_data, lat_data, res_factor=res_factor)
+
+ return _wrapper
+
+
+def _call_map_blocks_interp(func, lon_dask, lat_dask, res_factor):
+ new_row_chunks = tuple(x * res_factor for x in lon_dask.chunks[0])
+ new_col_chunks = tuple(x * res_factor for x in lon_dask.chunks[1])
+ wrapped_func = _map_blocks_handler(func)
+ res = da.map_blocks(wrapped_func, lon_dask, lat_dask, res_factor,
+ new_axis=[0],
+ chunks=(2, new_row_chunks, new_col_chunks),
+ dtype=lon_dask.dtype,
+ meta=np.empty((2, 2, 2), dtype=lon_dask.dtype))
+ return res[0], res[1]
+
+
+def _rechunk_lonlat_if_needed(lon_data, lat_data):
+ # take current chunk size and get a relatively similar chunk size
+ row_chunks = lon_data.chunks[0]
+ col_chunks = lon_data.chunks[1]
+ num_rows = lon_data.shape[0]
+ num_cols = lon_data.shape[-1]
+ good_row_chunks = all(x % ROWS_PER_SCAN == 0 for x in row_chunks)
+ good_col_chunks = len(col_chunks) == 1 and col_chunks[0] != num_cols
+ lonlat_same_chunks = lon_data.chunks == lat_data.chunks
+ if num_rows % ROWS_PER_SCAN != 0:
+ raise ValueError("Input longitude/latitude data does not consist of "
+ "whole scans (10 rows per scan).")
+ if good_row_chunks and good_col_chunks and lonlat_same_chunks:
+ return lon_data, lat_data
+
+ new_row_chunks = (row_chunks[0] // ROWS_PER_SCAN) * ROWS_PER_SCAN
+ lon_data = lon_data.rechunk((new_row_chunks, -1))
+ lat_data = lat_data.rechunk((new_row_chunks, -1))
+ return lon_data, lat_data
+
+
+def _map_blocks_handler(func):
+ def _map_blocks_wrapper(lon_array, lat_array, res_factor):
+ lons, lats = func(lon_array, lat_array, res_factor=res_factor)
+ return np.concatenate((lons[np.newaxis], lats[np.newaxis]), axis=0)
+ return _map_blocks_wrapper
+
+
+ at scanline_mapblocks
+def interpolate_geolocation_cartesian(lon_array, lat_array, res_factor=4):
+ """Interpolate MODIS navigation from 1000m resolution to 250m.
+
+ Python rewrite of the IDL function ``MODIS_GEO_INTERP_250`` but converts to cartesian (X, Y, Z) coordinates
+ first to avoid problems with the anti-meridian/poles.
+
+ Arguments:
+ lon_array: Longitude data as a 2D numpy, dask, or xarray DataArray object.
+ The input data is expected to represent 1000m geolocation.
+ lat_array: Latitude data as a 2D numpy, dask, or xarray DataArray object.
+ The input data is expected to represent 1000m geolocation.
+ res_factor (int): Expansion factor for the function. Should be 2 for
+ 500m output or 4 for 250m output.
+
+ Returns:
+ A two-element tuple (lon, lat).
+
+ """
+ num_rows, num_cols = lon_array.shape
+ num_scans = int(num_rows / ROWS_PER_SCAN)
+ x_in, y_in, z_in = lonlat2xyz(lon_array, lat_array)
+
+ # Create an array of indexes that we want our result to have
+ x = np.arange(res_factor * num_cols, dtype=np.float32) * (1. / res_factor)
+ # 0.375 for 250m, 0.25 for 500m
+ y = np.arange(res_factor * ROWS_PER_SCAN, dtype=np.float32) * \
+ (1. / res_factor) - (res_factor * (1. / 16) + (1. / 8))
+ x, y = np.meshgrid(x, y)
+ coordinates = np.array([y, x]) # Used by map_coordinates, major optimization
+
+ new_x = np.empty((num_rows * res_factor, num_cols * res_factor), dtype=lon_array.dtype)
+ new_y = new_x.copy()
+ new_z = new_x.copy()
+ nav_arrays = [(x_in, new_x), (y_in, new_y), (z_in, new_z)]
+
+ # Interpolate each scan, one at a time, otherwise the math doesn't work well
+ for scan_idx in range(num_scans):
+ # Calculate indexes
+ j0 = ROWS_PER_SCAN * scan_idx
+ j1 = j0 + ROWS_PER_SCAN
+ k0 = ROWS_PER_SCAN * res_factor * scan_idx
+ k1 = k0 + ROWS_PER_SCAN * res_factor
+
+ for nav_array, result_array in nav_arrays:
+ # Use bilinear interpolation for all 250 meter pixels
+ map_coordinates(nav_array[j0:j1, :], coordinates, output=result_array[k0:k1, :], order=1, mode='nearest')
+
+ if res_factor == 4:
+ # Use linear extrapolation for the first two 250 meter pixels along track
+ m, b = _calc_slope_offset_250(result_array, y, k0, 2)
+ result_array[k0 + 0, :] = m * y[0, 0] + b
+ result_array[k0 + 1, :] = m * y[1, 0] + b
+
+ # Use linear extrapolation for the last two 250 meter pixels along track
+ # m = (result_array[k0 + 37, :] - result_array[k0 + 34, :]) / (y[37, 0] - y[34, 0])
+ # b = result_array[k0 + 37, :] - m * y[37, 0]
+ m, b = _calc_slope_offset_250(result_array, y, k0, 34)
+ result_array[k0 + 38, :] = m * y[38, 0] + b
+ result_array[k0 + 39, :] = m * y[39, 0] + b
+ else:
+ # 500m
+ # Use linear extrapolation for the first two 250 meter pixels along track
+ m, b = _calc_slope_offset_500(result_array, y, k0, 1)
+ result_array[k0 + 0, :] = m * y[0, 0] + b
+
+ # Use linear extrapolation for the last two 250 meter pixels along track
+ m, b = _calc_slope_offset_500(result_array, y, k0, 17)
+ result_array[k0 + 19, :] = m * y[19, 0] + b
+
+ new_lons, new_lats = xyz2lonlat(new_x, new_y, new_z, low_lat_z=True)
+ return new_lons.astype(lon_array.dtype), new_lats.astype(lon_array.dtype)
+
+
+def _calc_slope_offset_250(result_array, y, start_idx, offset):
+ m = (result_array[start_idx + offset + 3, :] - result_array[start_idx + offset, :]) / \
+ (y[offset + 3, 0] - y[offset, 0])
+ b = result_array[start_idx + offset + 3, :] - m * y[offset + 3, 0]
+ return m, b
+
+
+def _calc_slope_offset_500(result_array, y, start_idx, offset):
+ m = (result_array[start_idx + offset + 1, :] - result_array[start_idx + offset, :]) / \
+ (y[offset + 1, 0] - y[offset, 0])
+ b = result_array[start_idx + offset + 1, :] - m * y[offset + 1, 0]
+ return m, b
+
+
+def modis_1km_to_250m(lon1, lat1):
+ """Interpolate MODIS geolocation from 1km to 250m resolution."""
+ return interpolate_geolocation_cartesian(lon1, lat1, res_factor=4)
+
+
+def modis_1km_to_500m(lon1, lat1):
+ """Interpolate MODIS geolocation from 1km to 500m resolution."""
+ return interpolate_geolocation_cartesian(lon1, lat1, res_factor=2)
=====================================
geotiepoints/tests/test_simple_modis_interpolator.py
=====================================
@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2021 Python-geotiepoints developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""Tests for simple MODIS interpolators."""
+
+import os
+import numpy as np
+import pytest
+import h5py
+import dask
+import dask.array as da
+import xarray as xr
+
+from geotiepoints.simple_modis_interpolator import modis_1km_to_250m, modis_1km_to_500m
+from .utils import CustomScheduler
+
+FILENAME_DATA = os.path.join(
+ os.path.dirname(__file__), '../../testdata/modis_test_data.h5')
+
+
+def _to_dask(arr):
+ return da.from_array(arr, chunks=4096)
+
+
+def _to_da(arr):
+ return xr.DataArray(_to_dask(arr), dims=['y', 'x'])
+
+
+def _load_h5_lonlat_vars(lon_var, lat_var):
+ h5f = h5py.File(FILENAME_DATA, 'r')
+ lon1 = h5f[lon_var]
+ lat1 = h5f[lat_var]
+ return lon1, lat1
+
+
+def _load_1km_lonlat_as_numpy():
+ lon1, lat1 = _load_h5_lonlat_vars('lon_1km', 'lat_1km')
+ return lon1[:], lat1[:]
+
+
+def _load_1km_lonlat_as_dask():
+ lon1, lat1 = _load_h5_lonlat_vars('lon_1km', 'lat_1km')
+ return _to_dask(lon1), _to_dask(lat1)
+
+
+def _load_1km_lonlat_as_xarray_dask():
+ lon1, lat1 = _load_h5_lonlat_vars('lon_1km', 'lat_1km')
+ return _to_da(lon1), _to_da(lat1)
+
+
+def _load_500m_lonlat_expected_as_xarray_dask():
+ h5f = h5py.File(FILENAME_DATA, 'r')
+ lon500 = _to_da(h5f['lon_500m'])
+ lat500 = _to_da(h5f['lat_500m'])
+ return lon500, lat500
+
+
+def _load_250m_lonlat_expected_as_xarray_dask():
+ h5f = h5py.File(FILENAME_DATA, 'r')
+ lon250 = _to_da(h5f['lon_250m'])
+ lat250 = _to_da(h5f['lat_250m'])
+ return lon250, lat250
+
+
+ at pytest.mark.parametrize(
+ ("input_func", "exp_func", "interp_func"),
+ [
+ (_load_1km_lonlat_as_xarray_dask, _load_500m_lonlat_expected_as_xarray_dask, modis_1km_to_500m),
+ (_load_1km_lonlat_as_xarray_dask, _load_250m_lonlat_expected_as_xarray_dask, modis_1km_to_250m),
+ (_load_1km_lonlat_as_dask, _load_500m_lonlat_expected_as_xarray_dask, modis_1km_to_500m),
+ (_load_1km_lonlat_as_dask, _load_250m_lonlat_expected_as_xarray_dask, modis_1km_to_250m),
+ (_load_1km_lonlat_as_numpy, _load_500m_lonlat_expected_as_xarray_dask, modis_1km_to_500m),
+ (_load_1km_lonlat_as_numpy, _load_250m_lonlat_expected_as_xarray_dask, modis_1km_to_250m),
+ ]
+)
+def test_basic_interp(input_func, exp_func, interp_func):
+ lon1, lat1 = input_func()
+ lons_exp, lats_exp = exp_func()
+
+ # when working with dask arrays, we shouldn't compute anything
+ with dask.config.set(scheduler=CustomScheduler(0)):
+ lons, lats = interp_func(lon1, lat1)
+
+ if hasattr(lons, "compute"):
+ lons, lats = da.compute(lons, lats)
+ # our "truth" values are from the modisinterpolator results
+ atol = 0.038 # 1e-2
+ rtol = 0
+ np.testing.assert_allclose(lons_exp, lons, atol=atol, rtol=rtol)
+ np.testing.assert_allclose(lats_exp, lats, atol=atol, rtol=rtol)
+ assert not np.any(np.isnan(lons))
+ assert not np.any(np.isnan(lats))
+
+
+def test_nonstandard_scan_size():
+ lon1, lat1 = _load_1km_lonlat_as_xarray_dask()
+ # remove 1 row from the end
+ lon1 = lon1[:-1]
+ lat1 = lat1[:-1]
+
+ pytest.raises(ValueError, modis_1km_to_250m, lon1, lat1)
=====================================
geotiepoints/tests/utils.py
=====================================
@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2021 Python-geotiepoints developers
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""Utilities for creating and checking tests."""
+
+
+class CustomScheduler(object):
+ """Scheduler raising an exception if data are computed too many times.
+
+ This only makes sense to include when working with dask-based arrays. To
+ use it::
+
+ with dask.config.set(scheduler=CustomScheduler(2)):
+ my_dask_arr.compute() # allowed
+ my_dask_arr.compute() # 2nd call, not allowed, fails
+
+ """
+
+ def __init__(self, max_computes=1):
+ """Set starting and maximum compute counts."""
+ self.max_computes = max_computes
+ self.total_computes = 0
+
+ def __call__(self, dsk, keys, **kwargs):
+ """Compute dask task and keep track of number of times we do so."""
+ import dask
+ self.total_computes += 1
+ if self.total_computes > self.max_computes:
+ raise RuntimeError("Too many dask computations were scheduled: "
+ "{}".format(self.total_computes))
+ return dask.get(dsk, keys, **kwargs)
=====================================
geotiepoints/version.py
=====================================
@@ -23,9 +23,9 @@ def get_keywords():
# setup.py/versioneer.py will grep for the variable names, so they must
# each be defined on a line of their own. _version.py will just call
# get_keywords().
- git_refnames = " (HEAD -> master, tag: v1.2.1)"
- git_full = "d083fad851b7dcbc0a397effff7b4623cc6099df"
- git_date = "2021-03-08 11:02:19 +0100"
+ git_refnames = " (tag: v1.3.0)"
+ git_full = "587d348b1711180b8e654ee650c08d4f6079df10"
+ git_date = "2021-09-12 14:01:30 -0500"
keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
return keywords
=====================================
setup.cfg
=====================================
@@ -4,6 +4,7 @@ release=1
[flake8]
max-line-length = 120
+ignore = D107
[versioneer]
VCS = git
=====================================
setup.py
=====================================
@@ -32,7 +32,7 @@ import numpy as np
from Cython.Build import cythonize
requirements = ['numpy', 'scipy', 'pandas']
-test_requires = ['h5py', 'xarray', 'dask']
+test_requires = ['pytest', 'pytest-cov', 'h5py', 'xarray', 'dask']
if sys.platform.startswith("win"):
extra_compile_args = []
@@ -71,7 +71,6 @@ if __name__ == "__main__":
cmdclass=cmdclass,
install_requires=requirements,
ext_modules=cythonize(EXTENSIONS),
- test_suite='geotiepoints.tests.suite',
tests_require=test_requires,
zip_safe=False
)
View it on GitLab: https://salsa.debian.org/debian-gis-team/python-geotiepoints/-/compare/21551727bf6e720f36672e8ba7828a990a78b12e...e5d72eb1b6745c9e45e0c47a1800c2e1930963f5
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/python-geotiepoints/-/compare/21551727bf6e720f36672e8ba7828a990a78b12e...e5d72eb1b6745c9e45e0c47a1800c2e1930963f5
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20210918/757ff05a/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list