Add cviirs-based fast modis interpolator

.travis.yml

@@ -4,7 +4,7 @@ env:
   - PYTHON_VERSION=$PYTHON_VERSION
   - NUMPY_VERSION=stable
   - MAIN_CMD='python setup.py'
-  - CONDA_DEPENDENCIES='cython pandas scipy numpy coveralls coverage h5py mock'
+  - CONDA_DEPENDENCIES='cython pandas scipy numpy coveralls coverage h5py mock dask xarray'
   - PIP_DEPENDENCIES=''
   - SETUP_XVFB=False
   - EVENT_TYPE='push pull_request'

appveyor.yml

@@ -3,7 +3,7 @@ environment:
     PYTHON: "C:\\conda"
     MINICONDA_VERSION: "latest"
     CMD_IN_ENV: "cmd /E:ON /V:ON /C .\\ci-helpers\\appveyor\\windows_sdk.cmd"
-    CONDA_DEPENDENCIES: "Cython pandas scipy numpy coveralls coverage h5py mock"
+    CONDA_DEPENDENCIES: "Cython pandas scipy numpy coveralls coverage h5py mock dask xarray"
     PIP_DEPENDENCIES: ""
     CONDA_CHANNELS: "conda-forge"
 

geotiepoints/modisinterpolator.py

@@ -0,0 +1,228 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2018 PyTroll community
+
+# Author(s):
+
+#   Martin Raspaud <[email protected]>
+
+# 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/>.
+
+"""Interpolation of geographical tiepoints using the second order interpolation
+scheme implemented in the CVIIRS software, as described here:
+Compact VIIRS SDR Product Format User Guide (V1J) 
+http://www.eumetsat.int/website/wcm/idc/idcplg?IdcService=GET_FILE&dDocName=PDF_DMT_708025&RevisionSelectionMethod=LatestReleased&Rendition=Web
+"""
+
+import xarray as xr
+import dask.array as da
+import numpy as np
+
+# TODO on interpolation:
+# - go over to cartesian coordinates for tricky situation (eg poles, dateline)
+
+R = 6371.
+# Aqua scan width and altitude in km
+scan_width = 10.00017
+H = 705.
+
+
+def compute_phi(zeta):
+    return np.arcsin(R * np.sin(zeta) / (R + H))
+
+
+def compute_theta(zeta, phi):
+    return zeta - phi
+
+
+def compute_zeta(phi):
+    return np.arcsin((R + H) * np.sin(phi) / R)
+
+
+def compute_expansion_alignment(satz_a, satz_b, satz_c, satz_d):
+    """All angles in radians."""
+    zeta_a = satz_a
+    zeta_b = satz_b
+
+    phi_a = compute_phi(zeta_a)
+    phi_b = compute_phi(zeta_b)
+    theta_a = compute_theta(zeta_a, phi_a)
+    theta_b = compute_theta(zeta_b, phi_b)
+    phi = (phi_a + phi_b) / 2
+    zeta = compute_zeta(phi)
+    theta = compute_theta(zeta, phi)
+
+    c_expansion = 4 * (((theta_a + theta_b) / 2 - theta) / (theta_a - theta_b))
+
+    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)
+
+    return c_expansion, c_alignment
+
+
+def get_corners(arr):
+    arr_a = arr[:, :-1, :-1]
+    arr_b = arr[:, :-1, 1:]
+    arr_c = arr[:, 1:, 1:]
+    arr_d = arr[:, 1:, :-1]
+    return arr_a, arr_b, arr_c, arr_d
+
+
+class ModisInterpolator():
+
+    def __init__(self, cres, fres):
+        if cres == 1000:
+            self.cscan_len = 10
+            self.cscan_width = 1
+            self.cscan_full_width = 1354
+        elif cres == 5000:
+            self.cscan_len = 2
+            self.cscan_width = 5
+            self.cscan_full_width = 271
+
+        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)
+        arr = da.concatenate((arr[:, :lines//2, :], arr, arr[:, -(lines//2):, :]), axis=1)
+        arr = da.repeat(arr.reshape((-1, self.cscan_full_width - 1)), cols, axis=1)
+        return da.hstack((arr, arr[:, -cols:]))
+
+    def _get_coords_1km(self, scans):
+        y = (np.arange((self.cscan_len + 1) * self.fscan_len) % self.fscan_len) + .5
+        y = y[self.fscan_len // 2:-(self.fscan_len // 2)]
+        y[:self.fscan_len//2] = np.arange(-self.fscan_len/2 + .5, 0)
+        y[-(self.fscan_len//2):] = np.arange(self.fscan_len + .5, self.fscan_len * 3 / 2)
+        y = np.tile(y, scans)
+
+        x = np.arange(self.fscan_full_width) % self.fscan_width
+        x[-self.fscan_width:] = np.arange(self.fscan_width, self.fscan_width * 2)
+        return x, y
+
+    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)
+        return da.hstack((arr[:, :2], arr, arr[:, -2:]))
+
+    def _get_coords_5km(self, scans):
+        y = np.arange(self.fscan_len * self.cscan_len) - 2
+        y = np.tile(y, scans)
+
+        x = (np.arange(self.fscan_full_width) - 2) % self.fscan_width
+        x[0] = -2
+        x[1] = -1
+        x[-2] = 5
+        x[-1] = 6
+        return x, y
+
+    def interpolate(self, lon1, lat1, satz1):
+        cscan_len = self.cscan_len
+        cscan_full_width = self.cscan_full_width
+
+        fscan_width = self.fscan_width
+        fscan_len = self.fscan_len
+
+        scans = lat1.shape[0] // cscan_len
+        latattrs = lat1.attrs
+        lonattrs = lon1.attrs
+        dims = lat1.dims
+        lat1 = lat1.data
+        lon1 = lon1.data
+        satz1 = satz1.data
+
+        lat1 = lat1.reshape((-1, cscan_len, cscan_full_width))
+        lon1 = lon1.reshape((-1, cscan_len, cscan_full_width))
+        satz1 = satz1.reshape((-1, cscan_len, cscan_full_width))
+
+        lats_a, lats_b, lats_c, lats_d = get_corners(lat1)
+        lons_a, lons_b, lons_c, lons_d = get_corners(lon1)
+        satz_a, satz_b, satz_c, satz_d = get_corners(da.deg2rad(satz1))
+        c_exp, c_ali = compute_expansion_alignment(satz_a, satz_b, satz_c, satz_d)
+
+        x, y = self.get_coords(scans)
+        i_rs, i_rt = da.meshgrid(x, y)
+
+        p_os = 0
+        p_ot = 0
+
+        s_s = (p_os + i_rs) * 1. / fscan_width
+        s_t = (p_ot + i_rt) * 1. / fscan_len
+
+        cols = fscan_width
+        lines = fscan_len
+
+        c_exp_full = self.expand_tiepoint_array(c_exp, lines, cols)
+        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)
+
+        lats_a = self.expand_tiepoint_array(lats_a, lines, cols)
+        lats_b = self.expand_tiepoint_array(lats_b, lines, cols)
+        lats_c = self.expand_tiepoint_array(lats_c, lines, cols)
+        lats_d = self.expand_tiepoint_array(lats_d, lines, cols)
+        lons_a = self.expand_tiepoint_array(lons_a, lines, cols)
+        lons_b = self.expand_tiepoint_array(lons_b, lines, cols)
+        lons_c = self.expand_tiepoint_array(lons_c, lines, cols)
+        lons_d = self.expand_tiepoint_array(lons_d, lines, cols)
+
+        lats_1 = (1 - a_scan) * lats_a + a_scan * lats_b
+        lats_2 = (1 - a_scan) * lats_d + a_scan * lats_c
+        lats = (1 - a_track) * lats_1 + a_track * lats_2
+
+        lons_1 = (1 - a_scan) * lons_a + a_scan * lons_b
+        lons_2 = (1 - a_scan) * lons_d + a_scan * lons_c
+        lons = (1 - a_track) * lons_1 + a_track * lons_2
+
+        return xr.DataArray(lons, attrs=lonattrs, dims=dims), xr.DataArray(lats, attrs=latattrs, dims=dims)
+
+
+def modis_1km_to_250m(lon1, lat1, satz1):
+
+    interp = ModisInterpolator(1000, 250)
+    return interp.interpolate(lon1, lat1, satz1)
+
+
+def modis_1km_to_500m(lon1, lat1, satz1):
+
+    interp = ModisInterpolator(1000, 500)
+    return interp.interpolate(lon1, lat1, satz1)
+
+
+def modis_5km_to_1km(lon1, lat1, satz1):
+
+    interp = ModisInterpolator(5000, 1000)
+    return interp.interpolate(lon1, lat1, satz1)

geotiepoints/tests/__init__.py

@@ -28,7 +28,8 @@
                                 test_modis,
                                 test_satelliteinterpolator,
                                 test_interpolator,
-                                test_multilinear)
+                                test_multilinear,
+                                test_modisinterpolator)
 
 
 import sys
@@ -58,6 +59,7 @@ def suite():
     mysuite.addTests(test_satelliteinterpolator.suite())
     mysuite.addTests(test_interpolator.suite())
     mysuite.addTests(test_multilinear.suite())
+    mysuite.addTests(test_modisinterpolator.suite())
     return mysuite
 
 if __name__ == '__main__':

geotiepoints/tests/test_modisinterpolator.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2018 Martin Raspaud
+
+# Author(s):
+
+#   Martin Raspaud <[email protected]>
+
+# 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/>.
+
+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
+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')
+        lon1 = to_da(h5f['lon_1km'])
+        lat1 = to_da(h5f['lat_1km'])
+        satz1 = to_da(h5f['satz_1km'])
+
+        lon250 = to_da(h5f['lon_250m'])
+        lon500 = to_da(h5f['lon_500m'])
+
+        lat250 = to_da(h5f['lat_250m'])
+        lat500 = to_da(h5f['lat_500m'])
+
+        lons, lats = modis_1km_to_250m(lon1, lat1, satz1)
+        self.assertTrue(np.allclose(lon250, lons, atol=1e-2))
+        self.assertTrue(np.allclose(lat250, lats, atol=1e-2))
+
+        lons, lats = modis_1km_to_500m(lon1, lat1, satz1)
+        self.assertTrue(np.allclose(lon500, lons, atol=1e-2))
+        self.assertTrue(np.allclose(lat500, lats, atol=1e-2))
+
+        lat5 = lat1[2::5, 2::5]
+        lon5 = lon1[2::5, 2::5]
+
+        satz5 = satz1[2::5, 2::5]
+        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))
+
+def suite():
+    """The suite for MODIS"""
+    loader = unittest.TestLoader()
+    mysuite = unittest.TestSuite()
+    mysuite.addTest(loader.loadTestsFromTestCase(TestModisInterpolator))
+
+    return mysuite
+
+if __name__ == "__main__":
+    unittest.main()

setup.py

@@ -38,7 +38,7 @@
 
 requirements = ['numpy', 'scipy', 'pandas'],
 # unittest2 is required by h5py 2.8.0rc:
-test_requires = ['h5py', 'unittest2']
+test_requires = ['h5py', 'unittest2', 'xarray', 'dask']
 
 if sys.platform.startswith("win"):
     extra_compile_args = []