[Git][debian-gis-team/python-geotiepoints][master] 8 commits: New upstream version 1.2.0
Antonio Valentino
gitlab at salsa.debian.org
Tue Jun 9 07:50:36 BST 2020
Antonio Valentino pushed to branch master at Debian GIS Project / python-geotiepoints
Commits:
7a9f813c by Antonio Valentino at 2020-06-09T05:35:17+00:00
New upstream version 1.2.0
- - - - -
9b8952c6 by Antonio Valentino at 2020-06-09T05:36:07+00:00
Update upstream source from tag 'upstream/1.2.0'
Update to upstream version '1.2.0'
with Debian dir bb89a1b284412b05e17e61efbb2e559c9476c738
- - - - -
1c551529 by Antonio Valentino at 2020-06-09T05:37:19+00:00
New upstrean release
- - - - -
eb68707c by Antonio Valentino at 2020-06-09T05:44:01+00:00
Update debian/copyright
- - - - -
d7f3f37c by Antonio Valentino at 2020-06-09T07:48:39+02:00
Referesh all patches
- - - - -
6337e3af by Antonio Valentino at 2020-06-09T06:13:19+00:00
Improve clean target
- - - - -
df9d2f4d by Antonio Valentino at 2020-06-09T06:43:21+00:00
Improve testing
- - - - -
f71482d0 by Antonio Valentino at 2020-06-09T06:43:58+00:00
Set distribution to unstable
- - - - -
16 changed files:
- + .stickler.yml
- .travis.yml
- + AUTHORS.md
- CHANGELOG.md
- appveyor.yml
- debian/changelog
- debian/copyright
- debian/patches/0001-Remove-dependency-form-unittest2.patch
- debian/patches/0002-Fix-tests.patch
- debian/rules
- geotiepoints/modisinterpolator.py
- geotiepoints/tests/test_modisinterpolator.py
- + geotiepoints/tests/test_viiinterpolator.py
- geotiepoints/version.py
- + geotiepoints/viiinterpolator.py
- setup.py
Changes:
=====================================
.stickler.yml
=====================================
@@ -0,0 +1,5 @@
+linters:
+ flake8:
+ python: 3
+ config: setup.cfg
+
=====================================
.travis.yml
=====================================
@@ -11,11 +11,6 @@ env:
- CONDA_CHANNELS='conda-forge'
matrix:
include:
- - env: PYTHON_VERSION=2.7
- os: linux
- - env: PYTHON_VERSION=2.7
- os: osx
- language: generic
- env: PYTHON_VERSION=3.6
os: linux
- env: PYTHON_VERSION=3.6
=====================================
AUTHORS.md
=====================================
@@ -0,0 +1,17 @@
+# Project Contributors
+
+The following people have made contributions to this project:
+
+<!--- Use your GitHub account or any other personal reference URL --->
+<!--- If you wish to not use your real name, please use your github username --->
+<!--- The list should be alphabetical by last name if possible, with github usernames at the bottom --->
+<!--- See https://gist.github.com/djhoese/52220272ec73b12eb8f4a29709be110d for auto-generating parts of this list --->
+
+- [Amit Aronovitch (AmitAronovitch)](https://github.com/AmitAronovitch)
+- [Adam Dybbroe (adybbroe)](https://github.com/adybbroe)
+- [Rolf Helge Pfeiffer (HelgeDMI)](https://github.com/HelgeDMI)
+- [David Hoese (djhoese)](https://github.com/djhoese)
+- [Mikhail Itkin (mitkin)](https://github.com/mitkin)
+- [Sauli Joro (sjoro)](https://github.com/sjoro)
+- [Panu Lahtinen (pnuu)](https://github.com/pnuu)
+- [Martin Raspaud (mraspaud)](https://github.com/mraspaud)
=====================================
CHANGELOG.md
=====================================
@@ -1,4 +1,21 @@
-## Version v1.1.8 (2019/04/24)
+## Version 1.2.0 (2020/06/05)
+
+
+### Pull Requests Merged
+
+#### Bugs fixed
+
+* [PR 19](https://github.com/pytroll/python-geotiepoints/pull/19) - Fix interpolation of symetrical tiepoints
+
+#### Features added
+
+* [PR 22](https://github.com/pytroll/python-geotiepoints/pull/22) - Add VII interpolator.
+* [PR 16](https://github.com/pytroll/python-geotiepoints/pull/16) - Add MODIS 5km to 500m and 250m interpolation
+
+In this release 3 pull requests were closed.
+
+
+## Version 1.1.8 (2019/04/24)
### Issues Closed
=====================================
appveyor.yml
=====================================
@@ -8,15 +8,14 @@ environment:
CONDA_CHANNELS: "conda-forge"
matrix:
- - PYTHON: "C:\\Python27_64"
- PYTHON_VERSION: "2.7"
- PYTHON_ARCH: "64"
- NUMPY_VERSION: "stable"
-
- PYTHON: "C:\\Python36_64"
PYTHON_VERSION: "3.6"
PYTHON_ARCH: "64"
NUMPY_VERSION: "stable"
+ - PYTHON: "C:\\Python37_64"
+ PYTHON_VERSION: "3.7"
+ PYTHON_ARCH: "64"
+ NUMPY_VERSION: "stable"
install:
- "git clone --depth 1 git://github.com/astropy/ci-helpers.git"
=====================================
debian/changelog
=====================================
@@ -1,10 +1,17 @@
-python-geotiepoints (1.1.8-3) UNRELEASED; urgency=medium
+python-geotiepoints (1.2.0-1) unstable; urgency=medium
+ * New upstream release.
* Bump Standards-Version to 4.5.0, no changes.
* debian/control:
- explicitly set Rules-Requires-Root: no
+ * Update debian/copyright.
+ * debian/patches:
+ - refresh all patches
+ - improve 0002-Fix-tests.patch
+ * debian/rules:
+ - improve clean target
- -- Antonio Valentino <antonio.valentino at tiscali.it> Mon, 30 Sep 2019 19:53:11 +0200
+ -- Antonio Valentino <antonio.valentino at tiscali.it> Tue, 09 Jun 2020 06:43:43 +0000
python-geotiepoints (1.1.8-2) unstable; urgency=medium
=====================================
debian/copyright
=====================================
@@ -4,11 +4,12 @@ Upstream-Contact: Adam Dybbroe, Martin Raspaud <martin.raspaud at smhi.se>
Source: https://github.com/pytroll/python-geotiepoints
Files: *
-Copyright: 2010-2018 PyTroll Team
+Copyright: 2010-2020 PyTroll Team
Martin Raspaud <martin.raspaud at smhi.se>
Adam Dybbroe <adam.dybbroe at smhi.se>
Panu Lahtinen <panu.lahtinen at fmi.fi>
PyTroll community
+ Python-geotiepoints developers
License: GPL-3+
Files: versioneer.py
@@ -16,7 +17,7 @@ Copyright: 2018, Brian Warner
License: CC0-1.0
Files: debian/*
-Copyright: 2018 Antonio Valentino <antonio.valentino at tiscali.it>
+Copyright: 2018-2020 Antonio Valentino <antonio.valentino at tiscali.it>
License: GPL-3+
License: GPL-3+
=====================================
debian/patches/0001-Remove-dependency-form-unittest2.patch
=====================================
@@ -8,7 +8,7 @@ It is not necessary for python3 packages
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/setup.py b/setup.py
-index d25f5e8..d6fa3f5 100644
+index 9a2eba4..61cc46d 100644
--- a/setup.py
+++ b/setup.py
@@ -33,7 +33,7 @@ import versioneer
=====================================
debian/patches/0002-Fix-tests.patch
=====================================
@@ -3,37 +3,31 @@ Date: Sun, 23 Dec 2018 10:36:27 +0000
Subject: Fix tests
---
- geotiepoints/tests/__init__.py | 23 ++++++++++++++++++-----
- setup.py | 2 +-
- 2 files changed, 19 insertions(+), 6 deletions(-)
+ geotiepoints/tests/__init__.py | 18 +++++++++++++++---
+ geotiepoints/tests/test_modis.py | 4 ++++
+ geotiepoints/tests/test_modisinterpolator.py | 1 +
+ setup.py | 2 +-
+ 4 files changed, 21 insertions(+), 4 deletions(-)
diff --git a/geotiepoints/tests/__init__.py b/geotiepoints/tests/__init__.py
-index 92150bc..b45394e 100644
+index 92150bc..63ea1e2 100644
--- a/geotiepoints/tests/__init__.py
+++ b/geotiepoints/tests/__init__.py
-@@ -25,11 +25,19 @@
- from geotiepoints import interpolator
-
- from geotiepoints.tests import (test_geointerpolator,
-- test_modis,
-+ # test_modis,
- test_satelliteinterpolator,
- test_interpolator,
+@@ -31,6 +31,13 @@ from geotiepoints.tests import (test_geointerpolator,
test_multilinear,
-- test_modisinterpolator)
-+ #test_modisinterpolator,
-+ )
-+
+ test_modisinterpolator)
+
+import os
+_DATADIR = os.path.join(os.path.dirname(__file__), '../../testdata')
+try:
+ import h5py
+except ImportError:
+ h5py = False
-
++
import sys
-@@ -49,17 +57,22 @@ def suite():
+ if sys.version_info < (2, 7):
+@@ -49,17 +56,22 @@ def suite():
mysuite = unittest.TestSuite()
if not TRAVIS:
# Test sphinx documentation pages:
@@ -59,8 +53,42 @@ index 92150bc..b45394e 100644
return mysuite
if __name__ == '__main__':
+diff --git a/geotiepoints/tests/test_modis.py b/geotiepoints/tests/test_modis.py
+index 2f921cb..464fbf9 100644
+--- a/geotiepoints/tests/test_modis.py
++++ b/geotiepoints/tests/test_modis.py
+@@ -47,6 +47,8 @@ class TestMODIS(unittest.TestCase):
+ def setUp(self):
+ pass
+
++ @unittest.skipIf(not os.path.isfile(FILENAME_FULL), 'deta file not available')
++ @unittest.skipIf(not os.path.isfile(FILENAME_5KM), 'deta file not available')
+ def test_5_to_1(self):
+ """test the 5km to 1km interpolation facility
+ """
+@@ -64,6 +66,8 @@ class TestMODIS(unittest.TestCase):
+ self.assert_(np.allclose(tlons, glons, atol=0.05))
+ self.assert_(np.allclose(tlats, glats, atol=0.05))
+
++ @unittest.skipIf(not os.path.isfile(FILENAME_250M_RESULT), 'deta file not available')
++ @unittest.skipIf(not os.path.isfile(FILENAME_250M_INPUT), 'deta file not available')
+ def test_1000m_to_250m(self):
+ """test the 1 km to 250 meter interpolation facility
+ """
+diff --git a/geotiepoints/tests/test_modisinterpolator.py b/geotiepoints/tests/test_modisinterpolator.py
+index f9f8bb5..ee295e0 100644
+--- a/geotiepoints/tests/test_modisinterpolator.py
++++ b/geotiepoints/tests/test_modisinterpolator.py
+@@ -40,6 +40,7 @@ def to_da(arr):
+ return xr.DataArray(da.from_array(arr, chunks=4096), dims=['y', 'x'])
+
+
++ at unittest.skipIf(not os.path.isfile(FILENAME_DATA), 'data file not available')
+ class TestModisInterpolator(unittest.TestCase):
+ def test_modis(self):
+ h5f = h5py.File(FILENAME_DATA, 'r')
diff --git a/setup.py b/setup.py
-index d6fa3f5..d7fa74a 100644
+index 61cc46d..8841c7b 100644
--- a/setup.py
+++ b/setup.py
@@ -128,7 +128,7 @@ if __name__ == "__main__":
@@ -71,4 +99,4 @@ index d6fa3f5..d7fa74a 100644
+ packages=['geotiepoints', 'geotiepoints.tests'],
# packages=find_packages(),
setup_requires=['numpy'],
- python_requires='>=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*',
+ python_requires='>=3.4',
=====================================
debian/rules
=====================================
@@ -16,6 +16,7 @@ override_dh_auto_clean:
dh_auto_clean
$(RM) geotiepoints/multilinear*.so
+ $(RM) -r python_geotiepoints.egg-info
override_dh_python3:
dh_python3
=====================================
geotiepoints/modisinterpolator.py
=====================================
@@ -29,6 +29,7 @@ http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_
import xarray as xr
import dask.array as da
import numpy as np
+import warnings
R = 6371.
# Aqua scan width and altitude in km
@@ -60,15 +61,17 @@ def compute_expansion_alignment(satz_a, satz_b, satz_c, satz_d):
phi = (phi_a + phi_b) / 2
zeta = compute_zeta(phi)
theta = compute_theta(zeta, phi)
+ # Workaround for tiepoints symetrical about the subsatellite-track
+ denominator = np.where(theta_a == theta_b, theta_a * 2, theta_a - theta_b)
- c_expansion = 4 * (((theta_a + theta_b) / 2 - theta) / (theta_a - theta_b))
+ c_expansion = 4 * (((theta_a + theta_b) / 2 - theta) / denominator)
sin_beta_2 = scan_width / (2 * H)
d = ((R + H) / R * np.cos(phi) - np.cos(zeta)) * sin_beta_2
e = np.cos(zeta) - np.sqrt(np.cos(zeta) ** 2 - d ** 2)
- c_alignment = 4 * e * np.sin(zeta) / (theta_a - theta_b)
+ c_alignment = 4 * e * np.sin(zeta) / denominator
return c_expansion, c_alignment
@@ -81,13 +84,15 @@ def get_corners(arr):
return arr_a, arr_b, arr_c, arr_d
-class ModisInterpolator():
+class ModisInterpolator(object):
def __init__(self, cres, fres, cscan_full_width=None):
if cres == 1000:
self.cscan_len = 10
self.cscan_width = 1
self.cscan_full_width = 1354
+ self.get_coords = self._get_coords_1km
+ self.expand_tiepoint_array = self._expand_tiepoint_array_1km
elif cres == 5000:
self.cscan_len = 2
self.cscan_width = 5
@@ -95,25 +100,21 @@ class ModisInterpolator():
self.cscan_full_width = 271
else:
self.cscan_full_width = cscan_full_width
+ self.expand_tiepoint_array = self._expand_tiepoint_array_5km
+ self.get_coords = self._get_coords_5km
if fres == 250:
self.fscan_width = 4 * self.cscan_width
self.fscan_full_width = 1354 * 4
self.fscan_len = 4 * 10 // self.cscan_len
- self.get_coords = self._get_coords_1km
- self.expand_tiepoint_array = self._expand_tiepoint_array_1km
elif fres == 500:
self.fscan_width = 2 * self.cscan_width
self.fscan_full_width = 1354 * 2
self.fscan_len = 2 * 10 // self.cscan_len
- self.get_coords = self._get_coords_1km
- self.expand_tiepoint_array = self._expand_tiepoint_array_1km
elif fres == 1000:
self.fscan_width = 1 * self.cscan_width
self.fscan_full_width = 1354
self.fscan_len = 1 * 10 // self.cscan_len
- self.get_coords = self._get_coords_5km
- self.expand_tiepoint_array = self._expand_tiepoint_array_5km
def _expand_tiepoint_array_1km(self, arr, lines, cols):
arr = da.repeat(arr, lines, axis=1)
@@ -135,10 +136,11 @@ class ModisInterpolator():
def _expand_tiepoint_array_5km(self, arr, lines, cols):
arr = da.repeat(arr, lines * 2, axis=1)
arr = da.repeat(arr.reshape((-1, self.cscan_full_width - 1)), cols, axis=1)
+ factor = self.fscan_width // self.cscan_width
if self.cscan_full_width == 271:
- return da.hstack((arr[:, :2], arr, arr[:, -2:]))
+ return da.hstack((arr[:, :2 * factor], arr, arr[:, -2 * factor:]))
else:
- return da.hstack((arr[:, :2], arr, arr[:, -5:], arr[:, -2:]))
+ return da.hstack((arr[:, :2 * factor], arr, arr[:, -self.fscan_width:], arr[:, -2 * factor:]))
def _get_coords_5km(self, scans):
y = np.arange(self.fscan_len * self.cscan_len) - 2
@@ -194,7 +196,7 @@ class ModisInterpolator():
c_ali_full = self.expand_tiepoint_array(c_ali, lines, cols)
a_track = s_t
- a_scan = (s_s + s_s * (1 - s_s) * c_exp_full + s_t*(1 - s_t) * c_ali_full)
+ a_scan = (s_s + s_s * (1 - s_s) * c_exp_full + s_t * (1 - s_t) * c_ali_full)
res = []
@@ -231,22 +233,39 @@ class ModisInterpolator():
def modis_1km_to_250m(lon1, lat1, satz1):
-
+ """Interpolate MODIS geolocation from 1km to 250m resolution."""
interp = ModisInterpolator(1000, 250)
return interp.interpolate(lon1, lat1, satz1)
def modis_1km_to_500m(lon1, lat1, satz1):
-
+ """Interpolate MODIS geolocation from 1km to 500m resolution."""
interp = ModisInterpolator(1000, 500)
return interp.interpolate(lon1, lat1, satz1)
def modis_5km_to_1km(lon1, lat1, satz1):
-
+ """Interpolate MODIS geolocation from 5km to 1km resolution."""
interp = ModisInterpolator(5000, 1000, lon1.shape[1])
return interp.interpolate(lon1, lat1, satz1)
+
+def modis_5km_to_500m(lon1, lat1, satz1):
+ """Interpolate MODIS geolocation from 5km to 500m resolution."""
+ warnings.warn("Interpolating 5km geolocation to 500m resolution "
+ "may result in poor quality")
+ interp = ModisInterpolator(5000, 500, lon1.shape[1])
+ return interp.interpolate(lon1, lat1, satz1)
+
+
+def modis_5km_to_250m(lon1, lat1, satz1):
+ """Interpolate MODIS geolocation from 5km to 250m resolution."""
+ warnings.warn("Interpolating 5km geolocation to 250m resolution "
+ "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
=====================================
geotiepoints/tests/test_modisinterpolator.py
=====================================
@@ -24,16 +24,22 @@ import unittest
import numpy as np
import h5py
import os
-from geotiepoints.modisinterpolator import modis_1km_to_250m, modis_1km_to_500m, modis_5km_to_1km
+from geotiepoints.modisinterpolator import (modis_1km_to_250m,
+ modis_1km_to_500m,
+ modis_5km_to_1km,
+ modis_5km_to_500m,
+ modis_5km_to_250m)
FILENAME_DATA = os.path.join(
os.path.dirname(__file__), '../../testdata/modis_test_data.h5')
+
def to_da(arr):
import xarray as xr
import dask.array as da
return xr.DataArray(da.from_array(arr, chunks=4096), dims=['y', 'x'])
+
class TestModisInterpolator(unittest.TestCase):
def test_modis(self):
h5f = h5py.File(FILENAME_DATA, 'r')
@@ -63,6 +69,20 @@ class TestModisInterpolator(unittest.TestCase):
self.assertTrue(np.allclose(lon1, lons, atol=1e-2))
self.assertTrue(np.allclose(lat1, lats, atol=1e-2))
+ # 5km to 500m
+ lons, lats = modis_5km_to_500m(lon5, lat5, satz5)
+ self.assertEqual(lon500.shape, lons.shape)
+ self.assertEqual(lat500.shape, lats.shape)
+ # self.assertTrue(np.allclose(lon500, lons, atol=1e-2))
+ # self.assertTrue(np.allclose(lat500, lats, atol=1e-2))
+
+ # 5km to 250m
+ lons, lats = modis_5km_to_250m(lon5, lat5, satz5)
+ self.assertEqual(lon250.shape, lons.shape)
+ self.assertEqual(lat250.shape, lats.shape)
+ # self.assertTrue(np.allclose(lon250, lons, atol=1e-2))
+ # self.assertTrue(np.allclose(lat250, lats, atol=1e-2))
+
# Test level 2
lat5 = lat1[2::5, 2:-5:5]
lon5 = lon1[2::5, 2:-5:5]
@@ -71,7 +91,13 @@ class TestModisInterpolator(unittest.TestCase):
lons, lats = modis_5km_to_1km(lon5, lat5, satz5)
self.assertTrue(np.allclose(lon1, lons, atol=1e-2))
self.assertTrue(np.allclose(lat1, lats, atol=1e-2))
-
+
+ # Test nans issue (#19)
+ satz1 = to_da(abs(np.linspace(-65.4, 65.4, 1354)).repeat(20).reshape(-1, 20).T)
+ lons, lats = modis_1km_to_500m(lon1, lat1, satz1)
+ self.assertFalse(np.any(np.isnan(lons.compute())))
+ self.assertFalse(np.any(np.isnan(lats.compute())))
+
def test_poles_datum(self):
import xarray as xr
h5f = h5py.File(FILENAME_DATA, 'r')
@@ -100,5 +126,6 @@ def suite():
return mysuite
+
if __name__ == "__main__":
unittest.main()
=====================================
geotiepoints/tests/test_viiinterpolator.py
=====================================
@@ -0,0 +1,252 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2017-2020 Python-geotiepoints developers
+#
+# This file is part of python-geotiepoints.
+#
+# python-geotiepoints 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.
+#
+# python-geotiepoints 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
+# python-geotiepoints. If not, see <http://www.gnu.org/licenses/>.
+
+"""Test of the interpolation of geographical tiepoints for the VII products.
+It follows the description provided in document "EPS-SG VII Level 1B Product Format Specification".
+"""
+
+import unittest
+import numpy as np
+import xarray as xr
+from geotiepoints.viiinterpolator import tie_points_interpolation, tie_points_geo_interpolation
+
+
+TEST_N_SCANS = 2
+TEST_TIE_POINTS_FACTOR = 2
+TEST_SCAN_ALT_TIE_POINTS = 3
+TEST_VALID_ALT_TIE_POINTS = TEST_SCAN_ALT_TIE_POINTS * TEST_N_SCANS
+TEST_INVALID_ALT_TIE_POINTS = TEST_SCAN_ALT_TIE_POINTS * TEST_N_SCANS + 1
+TEST_ACT_TIE_POINTS = 4
+
+# Results of latitude/longitude interpolation with simple interpolation on coordinates
+TEST_LON_1 = np.array(
+ [[-12., -11.5, -11., -10.5, -9., -8.5, -8., -7.5],
+ [-9., -8.5, -8., -7.5, -6., -5.5, -5., -4.5],
+ [-6., -5.5, -5., -4.5, -3., -2.5, -2., -1.5],
+ [-3., -2.5, -2., -1.5, 0., 0.5, 1., 1.5],
+ [0., 0.5, 1., 1.5, 3., 3.5, 4., 4.5],
+ [3., 3.5, 4., 4.5, 6., 6.5, 7., 7.5]]
+)
+TEST_LAT_1 = np.array(
+ [[0., 0.5, 1., 1.5, 3., 3.5, 4., 4.5],
+ [3., 3.5, 4., 4.5, 6., 6.5, 7., 7.5],
+ [6., 6.5, 7., 7.5, 9., 9.5, 10., 10.5],
+ [9., 9.5, 10., 10.5, 12., 12.5, 13., 13.5],
+ [12., 12.5, 13., 13.5, 15., 15.5, 16., 16.5],
+ [15., 15.5, 16., 16.5, 18., 18.5, 19., 19.5]]
+)
+
+# Results of latitude/longitude interpolation on cartesian coordinates (latitude above 60 degrees)
+TEST_LON_2 = np.array(
+ [[-12., -11.50003808, -11., -10.50011426, -9., -8.50026689, -8., -7.50034342],
+ [-9.00824726, -8.50989187, -8.01100418, -7.51272848, -6.01653996, -5.51842688, -5.01932226, -4.52129225],
+ [-6., -5.50049716, -5., -4.50057447, -3., -2.50073021, -2., -1.50080874],
+ [-3.02492451, -2.52706443, -2.02774808, -1.52997501, -0.03344942, 0.46414517, 0.96366893, 1.4611719],
+ [0., 0.49903263, 1., 1.49895241, 3., 3.49878988, 4., 4.49870746],
+ [2.9578336, 3.45514812, 3.9548757, 4.4520932, 5.94886832, 6.44588569, 6.94581415, 7.44272818]]
+)
+
+TEST_LAT_2 = np.array(
+ [[0., 0.49998096, 1., 1.49994287, 3., 3.49986656, 4., 4.4998283],
+ [2.99588485, 3.49506416, 3.99450923, 4.49364876, 5.99174708, 6.49080542, 6.99035883, 7.4893757],
+ [6., 6.49975143, 7., 7.49971278, 9., 9.49963492, 10., 10.49959566],
+ [8.98756357, 9.48649563, 9.98615477, 10.4850434, 11.98331018, 12.48210974, 12.98187246, 13.4806263],
+ [12., 12.49951634, 13., 13.49947623, 15., 15.49939498, 16., 16.49935377],
+ [14.97896116, 15.47762097, 15.97748548, 16.47609689, 17.97448854, 18.47300011, 18.97296495, 19.47142496]]
+)
+
+# Results of latitude/longitude interpolation on cartesian coordinates (longitude with a 360 degrees step)
+TEST_LON_3 = np.array(
+ [[-12., -11.50444038, -11., -10.50459822, -9., -8.50493209, -8., -7.50510905],
+ [-9.17477341, -8.68280267, -8.18102962, -7.68936161, -6.19433153, -5.70332248, -5.20141997, -4.71077058],
+ [-6., -5.50548573, -5., -4.50568668, -3., -2.50611746, -2., -1.506349],
+ [-3.2165963, -2.72673687, -2.22474246, -1.73531828, -0.24232275, 0.24613534, 0.7481615, 1.23608137],
+ [0., 0.49315061, 1., 1.49287934, 3., 3.49228746, 4., 4.49196335],
+ [2.72743411, 3.21414435, 3.71610443, 4.20213182, 5.69115252, 6.17562189, 6.67735289, 7.16092853]]
+)
+
+TEST_LAT_3 = np.array(
+ [[45., 45.49777998, 46., 46.49770107, 48., 48.49753416, 49., 49.49744569],
+ [47.91286617, 48.40886652, 48.90975264, 49.40560282, 50.90313463, 51.39865815, 51.8996091, 52.39495445],
+ [51., 51.49725738, 52., 52.49715691, 54., 54.50311196, 55., 55.50299846],
+ [54.11364234, 54.61463583, 55.10950687, 55.61043728, 57.10152529, 57.60232834, 58.09766771, 58.59840655],
+ [57., 57.50277937, 58., 58.50267347, 60., 60.50246826, 61., 61.5023687],
+ [60.09019289, 60.59080228, 61.0865661, 61.58711028, 63.07951189, 63.57992468, 64.07607638, 64.57642295]]
+)
+
+
+class TestViiInterpolator(unittest.TestCase):
+ """Test the vii_utils module."""
+
+ def setUp(self):
+ """Set up the test."""
+ # Create the arrays for the interpolation test
+ # The first has a valid number of n_tie_alt points (multiple of SCAN_ALT_TIE_POINTS)
+ self.valid_data_for_interpolation = xr.DataArray(
+ np.arange(
+ TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
+ dtype=np.float64,
+ ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
+ dims=('num_tie_points_act', 'num_tie_points_alt'),
+ )
+ # The second has an invalid number of n_tie_alt points (not multiple of SCAN_ALT_TIE_POINTS)
+ self.invalid_data_for_interpolation = xr.DataArray(
+ np.arange(
+ TEST_INVALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
+ dtype=np.float64,
+ ).reshape(TEST_ACT_TIE_POINTS, TEST_INVALID_ALT_TIE_POINTS),
+ dims=('num_tie_points_act', 'num_tie_points_alt'),
+ )
+ # Then two arrays containing valid longitude and latitude data
+ self.longitude = xr.DataArray(
+ np.linspace(
+ -12,
+ 11,
+ num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
+ dtype=np.float64,
+ ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
+ dims=('num_tie_points_act', 'num_tie_points_alt'),
+ )
+ self.latitude = xr.DataArray(
+ np.linspace(
+ 0,
+ 23,
+ num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
+ dtype=np.float64,
+ ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
+ dims=('num_tie_points_act', 'num_tie_points_alt'),
+ )
+ # Then one containing latitude data above 60 degrees
+ self.latitude_over60 = xr.DataArray(
+ np.linspace(
+ 45,
+ 68,
+ num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
+ dtype=np.float64,
+ ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS),
+ dims=('num_tie_points_act', 'num_tie_points_alt'),
+ )
+ # Then one containing longitude data with a 360 degrees step
+ self.longitude_over360 = xr.DataArray(
+ np.linspace(
+ -12,
+ 11,
+ num=TEST_VALID_ALT_TIE_POINTS * TEST_ACT_TIE_POINTS,
+ dtype=np.float64,
+ ).reshape(TEST_ACT_TIE_POINTS, TEST_VALID_ALT_TIE_POINTS) % 360.,
+ dims=('num_tie_points_act', 'num_tie_points_alt'),
+ )
+
+ def tearDown(self):
+ """Tear down the test."""
+ # Nothing to do
+ pass
+
+ def test_tie_points_interpolation(self):
+ """# Test the interpolation routine with valid and invalid input."""
+ # Test the interpolation routine with valid input
+ result_valid = tie_points_interpolation(
+ [self.valid_data_for_interpolation],
+ TEST_SCAN_ALT_TIE_POINTS,
+ TEST_TIE_POINTS_FACTOR
+ )[0]
+
+ act_points_interp = (TEST_ACT_TIE_POINTS - 1) * TEST_TIE_POINTS_FACTOR
+ num_scans = TEST_VALID_ALT_TIE_POINTS // TEST_SCAN_ALT_TIE_POINTS
+ scan_alt_points_interp = (TEST_SCAN_ALT_TIE_POINTS - 1) * TEST_TIE_POINTS_FACTOR
+
+ # It is easier to check the delta between interpolated points, which must be 1/8 of the original delta
+ # Across the track, it is possible to check the delta on the entire array
+ delta_axis_0 = 1.0 * TEST_VALID_ALT_TIE_POINTS / TEST_TIE_POINTS_FACTOR
+ self.assertTrue(np.allclose(
+ np.diff(result_valid, axis=0),
+ np.ones((act_points_interp - 1, num_scans * scan_alt_points_interp)) * delta_axis_0
+ ))
+
+ delta_axis_1 = 1.0 / TEST_TIE_POINTS_FACTOR
+ # Along the track, it is necessary to check the delta on each scan separately
+ for i in range(num_scans):
+ first_index = i*(TEST_SCAN_ALT_TIE_POINTS-1)*TEST_TIE_POINTS_FACTOR
+ last_index = (i+1)*(TEST_SCAN_ALT_TIE_POINTS-1)*TEST_TIE_POINTS_FACTOR
+ result_per_scan = result_valid[:, first_index:last_index]
+ self.assertTrue(np.allclose(
+ np.diff(result_per_scan, axis=1),
+ np.ones((act_points_interp, (TEST_SCAN_ALT_TIE_POINTS-1)*TEST_TIE_POINTS_FACTOR - 1)) * delta_axis_1
+ ))
+
+ self.assertEqual(len(result_valid.coords), 0)
+
+ # Test the interpolation routine with invalid input
+ with self.assertRaises(ValueError):
+ tie_points_interpolation(
+ [self.invalid_data_for_interpolation],
+ TEST_SCAN_ALT_TIE_POINTS,
+ TEST_TIE_POINTS_FACTOR
+ )[0]
+
+ def test_tie_points_geo_interpolation(self):
+ """# Test the coordinates interpolation routine with valid and invalid input."""
+ # Test the interpolation routine with valid input
+ lon, lat = tie_points_geo_interpolation(
+ self.longitude,
+ self.latitude,
+ TEST_SCAN_ALT_TIE_POINTS,
+ TEST_TIE_POINTS_FACTOR
+ )
+ self.assertTrue(np.allclose(lon, TEST_LON_1))
+ self.assertTrue(np.allclose(lat, TEST_LAT_1))
+
+ lon, lat = tie_points_geo_interpolation(
+ self.longitude_over360,
+ self.latitude,
+ TEST_SCAN_ALT_TIE_POINTS,
+ TEST_TIE_POINTS_FACTOR
+ )
+ self.assertTrue(np.allclose(lon, TEST_LON_2))
+ self.assertTrue(np.allclose(lat, TEST_LAT_2))
+
+ lon, lat = tie_points_geo_interpolation(
+ self.longitude,
+ self.latitude_over60,
+ TEST_SCAN_ALT_TIE_POINTS,
+ TEST_TIE_POINTS_FACTOR
+ )
+ self.assertTrue(np.allclose(lon, TEST_LON_3))
+ self.assertTrue(np.allclose(lat, TEST_LAT_3))
+
+ # Test the interpolation routine with invalid input (different dimensions of the two arrays)
+ with self.assertRaises(ValueError):
+ tie_points_geo_interpolation(
+ self.longitude,
+ self.invalid_data_for_interpolation,
+ TEST_SCAN_ALT_TIE_POINTS,
+ TEST_TIE_POINTS_FACTOR
+ )
+
+
+def suite():
+ """The suite for VII interpolator"""
+ loader = unittest.TestLoader()
+ mysuite = unittest.TestSuite()
+ mysuite.addTest(loader.loadTestsFromTestCase(TestViiInterpolator))
+
+ return mysuite
+
+
+if __name__ == "__main__":
+ unittest.main()
=====================================
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 = " (tag: v1.1.8)"
- git_full = "7c5cc8a887f8534cc2839c716c2c560aeaf77659"
- git_date = "2019-04-24 19:57:25 +0200"
+ git_refnames = " (HEAD -> master, tag: v1.2.0)"
+ git_full = "f28660ddfc898e07dc65bf5a031a5e8434aa3791"
+ git_date = "2020-06-05 12:34:19 +0200"
keywords = {"refnames": git_refnames, "full": git_full, "date": git_date}
return keywords
=====================================
geotiepoints/viiinterpolator.py
=====================================
@@ -0,0 +1,195 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Copyright (c) 2017-2020 Python-geotiepoints developers
+#
+# This file is part of python-geotiepoints.
+#
+# python-geotiepoints 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.
+#
+# python-geotiepoints 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
+# python-geotiepoints. If not, see <http://www.gnu.org/licenses/>.
+
+"""Interpolation of geographical tiepoints for the VII products.
+It follows the description provided in document "EPS-SG VII Level 1B Product Format Specification".
+Tiepoints are typically subsampled by a factor 8 with respect to the pixels, along and across the satellite track.
+Because of the bowtie effect, tiepoints at the scan edges are not continuous between neighbouring scans,
+therefore each scan has its own edge tiepoints in the along-track direction.
+Each scan typically extends on 3 tiepoints in the along-track direction.
+At the edges of a given scan (both along and across track) the tie points lie outside the original data point raster
+and are therefore excluded from the interpolation grid.
+The interpolation functions are implemented for xarray.DataArrays as input.
+"""
+
+import xarray as xr
+import dask.array as da
+import numpy as np
+
+# MEAN EARTH RADIUS AS DEFINED BY IUGG
+MEAN_EARTH_RADIUS = 6371008.7714 # [m]
+
+
+def tie_points_interpolation(data_on_tie_points, scan_alt_tie_points, tie_points_factor):
+ """Interpolate the data from the tie points to the pixel points.
+
+ The data are provided as a list of xarray DataArray objects, allowing to interpolate on several arrays
+ at the same time; however the individual arrays must have exactly the same dimensions.
+
+ Args:
+ data_on_tie_points: list of xarray DataArray objects containing the values defined on the tie points.
+ scan_alt_tie_points: number of tie points along the satellite track for each scan
+ tie_points_factor: sub-sampling factor of tie points wrt pixel points
+
+ Returns:
+ list of xarray DataArray objects containing the interpolated values on the pixel points.
+
+ """
+ # Extract the dimensions of the tie points array across and along track
+ n_tie_act, n_tie_alt = data_on_tie_points[0].shape
+ dim_act, dim_alt = data_on_tie_points[0].dims
+
+ # Check that the number of tie points along track is multiple of the number of tie points per scan
+ if n_tie_alt % scan_alt_tie_points != 0:
+ raise ValueError("The number of tie points in the along-route dimension must be a multiple of %d",
+ scan_alt_tie_points)
+
+ # Compute the number of scans
+ n_scans = n_tie_alt // scan_alt_tie_points
+
+ # Compute the dimensions of the pixel points array across and along track
+ n_pixel_act = (n_tie_act - 1) * tie_points_factor
+ n_pixel_alt = (n_tie_alt - 1) * tie_points_factor
+
+ # Create the grids used for interpolation across the track
+ tie_grid_act = da.arange(0, n_pixel_act + 1, tie_points_factor)
+ pixels_grid_act = da.arange(0, n_pixel_act)
+
+ # Create the grids used for the interpolation along the track (must not include the spurious points between scans)
+ tie_grid_alt = da.arange(0, n_pixel_alt + 1, tie_points_factor)
+ n_pixel_alt_per_scan = (scan_alt_tie_points - 1) * tie_points_factor
+ pixel_grid_alt = []
+ for j_scan in range(n_scans):
+ start_index_scan = j_scan * scan_alt_tie_points * tie_points_factor
+ pixel_grid_alt.append(da.arange(start_index_scan, start_index_scan + n_pixel_alt_per_scan))
+ pixel_grid_alt = da.concatenate(pixel_grid_alt)
+
+ # Loop on all arrays
+ data_on_pixel_points = []
+ for data in data_on_tie_points:
+
+ if data.shape != (n_tie_act, n_tie_alt) or data.dims != (dim_act, dim_alt):
+ raise ValueError("The dimensions of the arrays are not consistent")
+
+ # Interpolate using the xarray interp function twice: first across, then along the scan
+ # (much faster than interpolating directly in the two dimensions)
+ data = data.assign_coords({dim_act: tie_grid_act, dim_alt: tie_grid_alt})
+ data_pixel = data.interp({dim_act: pixels_grid_act}, assume_sorted=True) \
+ .interp({dim_alt: pixel_grid_alt}, assume_sorted=True).drop_vars([dim_act, dim_alt])
+
+ data_on_pixel_points.append(data_pixel)
+
+ return data_on_pixel_points
+
+
+def tie_points_geo_interpolation(longitude, latitude,
+ scan_alt_tie_points, tie_points_factor,
+ lat_threshold_use_cartesian=60.,
+ z_threshold_use_xy=0.8):
+ """Interpolate the geographical position from the tie points to the pixel points.
+
+ The longitude and latitude values are provided as xarray DataArray objects.
+
+ Args:
+ data_on_tie_points: list of xarray DataArray objects containing the values defined on the tie points.
+ scan_alt_tie_points: number of tie points along the satellite track for each scan
+ tie_points_factor: sub-sampling factor of tie points wrt pixel points
+
+ Returns:
+ list of xarray DataArray objects containing the interpolated values on the pixel points.
+
+ Args:
+ longitude: xarray DataArray containing the longitude values defined on the tie points (degrees).
+ latitude: xarray DataArray containing the latitude values defined on the tie points (degrees).
+ scan_alt_tie_points: number of tie points along the satellite track for each scan.
+ tie_points_factor: sub-sampling factor of tie points wrt pixel points.
+ lat_threshold_use_cartesian: latitude threshold to use cartesian coordinates.
+ z_threshold_use_xy: z threshold to compute latitude from x and y in cartesian coordinates.
+
+ Returns:
+ two xarray DataArray objects containing the interpolated longitude and latitude values on the pixel points.
+
+ """
+ # Check that the two arrays have the same dimensions
+ if longitude.shape != latitude.shape:
+ raise ValueError("The dimensions of longitude and latitude don't match")
+
+ # Determine if the interpolation should be done in cartesian or geodetic coordinates
+ to_cart = np.max(np.fabs(latitude)) > lat_threshold_use_cartesian or (np.max(longitude) - np.min(longitude)) > 180.
+
+ if to_cart:
+
+ x, y, z = _lonlat2xyz(longitude, latitude)
+
+ interp_x, interp_y, interp_z = tie_points_interpolation([x, y, z],
+ scan_alt_tie_points,
+ tie_points_factor)
+
+ interp_longitude, interp_latitude = _xyz2lonlat(interp_x, interp_y, interp_z, z_threshold_use_xy)
+
+ else:
+
+ interp_longitude, interp_latitude = tie_points_interpolation([longitude, latitude],
+ scan_alt_tie_points,
+ tie_points_factor)
+
+ return interp_longitude, interp_latitude
+
+
+def _lonlat2xyz(lons, lats):
+ """Convert longitudes and latitudes to cartesian coordinates.
+
+ Args:
+ lons: array containing the longitude values in degrees.
+ lats: array containing the latitude values in degrees.
+
+ Returns:
+ tuple of arrays containing the x, y, and z values in meters.
+
+ """
+ lons_rad = np.deg2rad(lons)
+ lats_rad = np.deg2rad(lats)
+ x_coords = MEAN_EARTH_RADIUS * np.cos(lats_rad) * np.cos(lons_rad)
+ y_coords = MEAN_EARTH_RADIUS * np.cos(lats_rad) * np.sin(lons_rad)
+ z_coords = MEAN_EARTH_RADIUS * np.sin(lats_rad)
+ return x_coords, y_coords, z_coords
+
+
+def _xyz2lonlat(x_coords, y_coords, z_coords, z_threshold_use_xy=0.8):
+ """Get longitudes and latitudes from cartesian coordinates.
+
+ Args:
+ x_coords: array containing the x values in meters.
+ y_coords: array containing the y values in meters.
+ z_coords: array containing the z values in meters.
+ z_threshold_use_xy: z threshold to compute latitude from x and y in cartesian coordinates.
+
+ Returns:
+ tuple of arrays containing the longitude and latitude values in degrees.
+
+ """
+ r = np.sqrt(x_coords ** 2 + y_coords ** 2)
+ thr_z = z_threshold_use_xy * MEAN_EARTH_RADIUS
+ lons = np.rad2deg(np.arccos(x_coords / r)) * np.sign(y_coords)
+ # Compute latitude from z at low z and from x and y at high z
+ lats = xr.where(
+ np.logical_and(np.less(z_coords, thr_z), np.greater(z_coords, -thr_z)),
+ 90. - np.rad2deg(np.arccos(z_coords / MEAN_EARTH_RADIUS)),
+ np.sign(z_coords) * (90. - np.rad2deg(np.arcsin(r / MEAN_EARTH_RADIUS)))
+ )
+ return lons, lats
=====================================
setup.py
=====================================
@@ -131,7 +131,7 @@ if __name__ == "__main__":
packages=['geotiepoints'],
# packages=find_packages(),
setup_requires=['numpy'],
- python_requires='>=2.7,!=3.0.*,!=3.1.*,!=3.2.*,!=3.3.*',
+ python_requires='>=3.4',
cmdclass=cmdclass,
install_requires=requirements,
View it on GitLab: https://salsa.debian.org/debian-gis-team/python-geotiepoints/-/compare/9a56c41e2624d2738fb014c398e91b02a37447c7...f71482d0629f0ae28878f418755cba5534705d62
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/python-geotiepoints/-/compare/9a56c41e2624d2738fb014c398e91b02a37447c7...f71482d0629f0ae28878f418755cba5534705d62
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/20200609/06e467d5/attachment-0001.html>
More information about the Pkg-grass-devel
mailing list