Added cartesian modis 1km to 250m interpolation method

py/polar2grid_modis/polar2grid/modis/modis_geo_interp_250.py

@@ -37,7 +37,7 @@
 :license:      GNU GPLv3
 """
 
-import numpy
+import numpy as np
 from scipy.ndimage.interpolation import map_coordinates
 
 import logging
@@ -48,6 +48,91 @@
 ROWS_PER_SCAN = 10
 # If we are going from 1000m to 250m we have 4 times the size of the original
 RES_FACTOR = 4
+EARTH_RADIUS = 6370997.0
+
+# FUTURE: Add this to the pytroll python-geotiepoints package if it can be more generalized
+# Most of the cartesian conversions were taken from the existing python-geotiepoints develop branch
+
+def get_lons_from_cartesian(x__, y__):
+    """Get longitudes from cartesian coordinates.
+    """
+    return np.rad2deg(np.arccos(x__ / np.sqrt(x__ ** 2 + y__ ** 2))) * np.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 - np.rad2deg(np.arccos(z__/EARTH_RADIUS)),
+                    np.sign(z__) *
+                    (90 - np.rad2deg(np.arcsin(np.sqrt(x__ ** 2 + y__ ** 2)
+                                         / EARTH_RADIUS))))
+    return lats
+
+
+def interpolate_geolocation_cartesian(lon_array, lat_array):
+    """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.
+
+    :param nav_array: MODIS 1km latitude array or 1km longitude array
+
+    :returns: MODIS 250m latitude array or 250m longitude array
+    """
+    num_rows,num_cols = lon_array.shape
+    num_scans = num_rows / ROWS_PER_SCAN
+
+    lons_rad = np.radians(lon_array)
+    lats_rad = np.radians(lat_array)
+    x_in = EARTH_RADIUS * np.cos(lats_rad) * np.cos(lons_rad)
+    y_in = EARTH_RADIUS * np.cos(lats_rad) * np.sin(lons_rad)
+    z_in = EARTH_RADIUS * np.sin(lats_rad)
+
+    # Create an array of indexes that we want our result to have
+    x = np.arange(RES_FACTOR * num_cols, dtype=np.float32) * 0.25
+    y = np.arange(RES_FACTOR * ROWS_PER_SCAN, dtype=np.float32) * 0.25 - 0.375
+    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=np.float32 )
+    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')
+
+            # Use linear extrapolation for the first two 250 meter pixels along track
+            m = (result_array[ k0 + 5, : ] - result_array[ k0 + 2, : ]) / (y[5,0] - y[2,0])
+            b = result_array[ k0 + 5, : ] - m * y[5,0]
+            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]
+            result_array[ k0 + 38, : ] = m * y[38,0] + b
+            result_array[ k0 + 39, : ] = m * y[39,0] + b
+
+    # Convert from cartesian to lat/lon space
+    new_lons = get_lons_from_cartesian(new_x, new_y)
+    new_lats = get_lats_from_cartesian(new_x, new_y, new_z)
+
+    return new_lons, new_lats
 
 
 def interpolate_geolocation(nav_array):
@@ -63,11 +148,11 @@ def interpolate_geolocation(nav_array):
     num_scans = num_rows / ROWS_PER_SCAN
 
     # Make a resulting array that is the right size for the new resolution
-    result_array = numpy.empty( (num_rows * RES_FACTOR, num_cols * RES_FACTOR), dtype=numpy.float32 )
-    x = numpy.arange(RES_FACTOR * num_cols, dtype=numpy.float32) * 0.25
-    y = numpy.arange(RES_FACTOR * ROWS_PER_SCAN, dtype=numpy.float32) * 0.25 - 0.375
-    x,y = numpy.meshgrid(x,y)
-    coordinates = numpy.array([y,x]) # Used by map_coordinates, major optimization
+    result_array = np.empty( (num_rows * RES_FACTOR, num_cols * RES_FACTOR), dtype=np.float32 )
+    x = np.arange(RES_FACTOR * num_cols, dtype=np.float32) * 0.25
+    y = np.arange(RES_FACTOR * ROWS_PER_SCAN, dtype=np.float32) * 0.25 - 0.375
+    x,y = np.meshgrid(x,y)
+    coordinates = np.array([y,x]) # Used by map_coordinates, major optimization
 
     # Interpolate each scan, one at a time, otherwise the math doesn't work well
     for scan_idx in range(num_scans):

py/polar2grid_modis/polar2grid/modis/modis_guidebook.py

@@ -40,8 +40,7 @@
 __docformat__ = "restructuredtext en"
 
 from polar2grid.core.frontend_utils import BaseFileReader, BaseMultiFileReader
-#from polar2grid.modis.modis_geo_interp_250 import interpolate_geolocation
-from geotiepoints import modis1kmto250m
+from polar2grid.modis.modis_geo_interp_250 import interpolate_geolocation_cartesian
 
 import os
 import logging
@@ -472,7 +471,7 @@ def get_swath_data(self, item, fill=None):
             if mask is not None:
                 data[mask] = numpy.nan
 
-            if None not in self.nav_interpolation["250"]:
+            if self.nav_interpolation["250"][0] is not None and self.nav_interpolation["250"][1] is not None:
                 LOG.debug("Returning previously interpolated 250m resolution geolocation data")
                 data = self.nav_interpolation["250"][not (item == K_LONGITUDE_250)]
                 self.nav_interpolation["250"] = [None, None]
@@ -483,7 +482,7 @@ def get_swath_data(self, item, fill=None):
             else:
                 self.nav_interpolation["250"][1] = data
 
-            if None in self.nav_interpolation["250"]:
+            if self.nav_interpolation["250"][0] is None or self.nav_interpolation["250"][1] is None:
                 # We don't have the other coordinate data yet
                 self.get_swath_data(K_LONGITUDE_250 if item == K_LATITUDE_250 else K_LATITUDE_250, fill=fill)
             else:
@@ -494,13 +493,7 @@ def get_swath_data(self, item, fill=None):
             LOG.info("Interpolating to higher resolution: %s" % (var_info.var_name,))
             lon_data, lat_data = self.nav_interpolation["250"]
 
-            if (lon_data.max() - lon_data.min()) > 180.0:
-                # we must be crossing the dateline or a pole, use the long running interpolation
-                LOG.info("Data crosses the dateline, interpolation will take longer than usual...")
-                new_lon_data, new_lat_data = modis1kmto250m(lon_data, lat_data)
-            else:
-                new_lon_data = interpolate_geolocation(lon_data)
-                new_lat_data = interpolate_geolocation(lat_data)
+            new_lon_data, new_lat_data = interpolate_geolocation_cartesian(lon_data, lat_data)
 
             new_lon_data[numpy.isnan(new_lon_data)] = fill
             new_lat_data[numpy.isnan(new_lat_data)] = fill