[Git][debian-gis-team/python-geotiepoints][upstream] New upstream version 1.2.0

Antonio Valentino gitlab at salsa.debian.org
Tue Jun 9 07:50:44 BST 2020



Antonio Valentino pushed to branch upstream 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
- - - - -


11 changed files:

- + .stickler.yml
- .travis.yml
- + AUTHORS.md
- CHANGELOG.md
- appveyor.yml
- 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"


=====================================
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/-/commit/7a9f813c04ab5923847341b5c9201092826fae85

-- 
View it on GitLab: https://salsa.debian.org/debian-gis-team/python-geotiepoints/-/commit/7a9f813c04ab5923847341b5c9201092826fae85
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/e8b40d81/attachment-0001.html>


More information about the Pkg-grass-devel mailing list