Merge pull request #143 from mhearne-usgs/interpfix

Added new global functions for converting geodicts to affine objects …

mapio/geodict.py

@@ -1,8 +1,10 @@
 #!/usr/bin/env python
 
+# third party imports
 import numpy as np
 from .dataset import DataSetException
 from osgeo import osr
+from rasterio.transform import Affine
 
 
 class GeoDict(object):
@@ -857,3 +859,41 @@ def validate(self, adjust):
             self._xmin += 360
         # if self._xmin == 180.0 and self._xmax < 0:
         #     self._xmin = -180.0
+
+
+def geodict_from_src(src):
+    """Return a geodict from a dataset object.
+
+    Args:
+        src (DatasetReader): Open rasterio DatasetReader object.
+    Returns:
+        GeoDict: GeoDict describing the DatasetReader object.
+    """
+    affine = src.transform
+    nx = src.width
+    ny = src.height
+    return geodict_from_affine(affine, ny, nx)
+
+
+def geodict_from_affine(affine, ny, nx):
+    geodict = {}
+    geodict['dx'] = affine.a
+    geodict['dy'] = -1 * affine.e
+    geodict['xmin'] = affine.xoff + geodict['dx'] / 2.0
+    geodict['ymax'] = affine.yoff - geodict['dy'] / 2.0
+    geodict['ny'] = ny
+    geodict['nx'] = nx
+    geodict['xmax'] = geodict['xmin'] + (geodict['nx'] - 1) * geodict['dx']
+    geodict['ymin'] = geodict['ymax'] - (geodict['ny'] - 1) * geodict['dy']
+
+    gd = GeoDict(geodict)
+    return gd
+
+
+def affine_from_geodict(geodict):
+    xoff = geodict.xmin - geodict.dx / 2.0
+    yoff = geodict.ymax + geodict.dy / 2.0
+    xres = geodict.dx
+    yres = -1 * geodict.dy
+    src = Affine.translation(xoff, yoff) * Affine.scale(xres, yres)
+    return src

mapio/grid2d.py

@@ -7,12 +7,11 @@
 import abc
 import textwrap
 import sys
-import re
 
 # third party imports
-from .gridbase import Grid
-from .dataset import DataSetException
-from .geodict import GeoDict
+from mapio.gridbase import Grid
+from mapio.dataset import DataSetException
+from mapio.geodict import GeoDict, affine_from_geodict
 
 import numpy as np
 from scipy import interpolate
@@ -1067,6 +1066,17 @@ def interpolate2(self, geodict, method='linear'):
           A new instance of the Grid2D class or subclass with interpolated
           data.
         """
+        if not self._geodict.contains(geodict):
+            msg = ('Input geodict not fully contained by host geodict.')
+            raise Exception(msg)
+
+        if self.getGeoDict().isAligned(geodict):
+            # this grid is already aligned, so let's just cut it at
+            # the input bounds
+            xmin, xmax, ymin, ymax = (geodict.xmin, geodict.xmax,
+                                      geodict.ymin, geodict.ymax)
+            newgrid = self.cut(xmin, xmax, ymin, ymax)
+            return newgrid
         resampling = None
         if method == 'linear':
             resampling = Resampling.bilinear
@@ -1077,31 +1087,26 @@ def interpolate2(self, geodict, method='linear'):
         else:
             raise DataSetException('Unknown interpolation method %s' % method)
 
-        destination = np.zeros((geodict.ny, geodict.nx))
-        src_transform = Affine.from_gdal(self._geodict.xmin -
-                                         self._geodict.dx / 2.0,
-                                         self._geodict.dx,
-                                         0.0,  # x rotation, not used by us
-                                         self._geodict.ymax +
-                                         self._geodict.dy / 2.0,
-                                         0.0,  # y rotation, not used by us
-                                         -1 * self._geodict.dy)  # their dy
-        # is negative
-        dst_transform = Affine.from_gdal(geodict.xmin - geodict.dx / 2.0,
-                                         geodict.dx,
-                                         0.0,  # x rotation, not used by us
-                                         geodict.ymax + geodict.dy / 2.0,
-                                         0.0,  # y rotation, not used by us
-                                         -1 * geodict.dy)  # their dy is negative
+        # we want the output data type to be 32 bit float in all cases
+        # except where the input is 64 bit float, in which case we need to
+        # preserve the precision.
+        dtype = np.float32
+        if self._data.dtype == np.float64:
+            dtype = np.float64
+
+        destination = np.zeros((geodict.ny, geodict.nx),
+                               dtype=dtype)
+        src_transform = affine_from_geodict(self._geodict)
+        dst_transform = affine_from_geodict(geodict)
         src_crs = CRS().from_string(GeoDict.DEFAULT_PROJ4).to_dict()
         dst_crs = CRS().from_string(GeoDict.DEFAULT_PROJ4).to_dict()
-        reproject(self._data.astype(np.float64), destination,
+        reproject(self._data.astype(dtype), destination,
                   src_transform=src_transform,
                   src_crs=src_crs,
                   src_nodata=np.nan,
                   dst_transform=dst_transform,
                   dst_crs=dst_crs,
-                  dst_nodata=np.nan,
+                  dst_nodata=float('nan'),
                   resampling=resampling)
 
         return self.__class__(destination, geodict)

mapio/reader.py

@@ -5,7 +5,7 @@
 
 # local imports
 from .grid2d import Grid2D
-from .geodict import GeoDict
+from .geodict import GeoDict, geodict_from_src
 
 
 def _is_hdf(filename):
@@ -31,29 +31,6 @@ def _is_hdf(filename):
     return is_hdf
 
 
-def _get_geodict_from_src(src):
-    """Return a geodict from a dataset object.
-
-    Args:
-        src (DatasetReader): Open rasterio DatasetReader object.
-    Returns:
-        GeoDict: GeoDict describing the DatasetReader object.
-    """
-    affine = src.transform
-    geodict = {}
-    geodict['dx'] = affine.a
-    geodict['dy'] = -1 * affine.e
-    geodict['xmin'] = affine.xoff + geodict['dx'] / 2.0
-    geodict['ymax'] = affine.yoff - geodict['dy'] / 2.0
-    geodict['ny'] = src.height
-    geodict['nx'] = src.width
-    geodict['xmax'] = geodict['xmin'] + (geodict['nx'] - 1) * geodict['dx']
-    geodict['ymin'] = geodict['ymax'] - (geodict['ny'] - 1) * geodict['dy']
-
-    gd = GeoDict(geodict)
-    return gd
-
-
 def _get_geodict_from_window(affine, window, data):
     """Return a geodict from a rasterio Window object.
 
@@ -229,6 +206,7 @@ def _read_data(src, samplegeodict, resample, method):
                                     pad=pad)
         data, src = _read_pixels(src, window)
         gd = _get_geodict_from_window(affine, window, data)
+
         gd.nodata = src.nodata
         grid = Grid2D(data, gd)
         return grid
@@ -255,6 +233,19 @@ def _read_data(src, samplegeodict, resample, method):
         left_block, src = _read_pixels(src, left_window)
 
         right_window = _geodict_to_window(sample_right, src, pad=pad)
+
+        # if the right window is a different width than expected,
+        # adjust the width to match the input number of columns
+        dwidth = (left_window.width + right_window.width) - samplegeodict.nx
+        if dwidth != 0:
+            col_off = right_window.col_off
+            row_off = right_window.row_off
+            width = right_window.width - dwidth
+            height = right_window.height
+            right_window = rasterio.windows.Window(col_off,
+                                                   row_off,
+                                                   width,
+                                                   height)
         right_block, src = _read_pixels(src, right_window)
 
         left_gd = _get_geodict_from_window(affine,
@@ -288,13 +279,13 @@ def get_file_geodict(filename):
         GeoDict: GeoDict: GeoDict describing the entire file.
     """
     src = rasterio.open(filename)
-    gd = _get_geodict_from_src(src)
+    gd = geodict_from_src(src)
     return gd
 
 
 def read(filename, samplegeodict=None, resample=False,
          method='linear', doPadding=False, padValue=np.nan,
-         apply_nan=True, force_cast=True):
+         apply_nan=True, force_cast=True, interp_approach='scipy'):
     """Read part or all of a rasterio file, resampling and padding as necessary.
 
     If samplegeodict is not provided, then the entire file will be read.
@@ -319,7 +310,8 @@ def read(filename, samplegeodict=None, resample=False,
                           from file.
         padValue (float): Value to insert in ring of padding pixels.
         apply_nan (bool): Convert nodata values to NaNs, upcasting to float if necessary.
-        force_cast (bool): If data values exceed 
+        force_cast (bool): If data values exceed range of values, cast upward.
+        interp_approach (str): One of 'scipy', 'rasterio'.
     """
     # use rasterio to read all formats
     src = rasterio.open(filename)
@@ -328,7 +320,7 @@ def read(filename, samplegeodict=None, resample=False,
     # if not, read the whole file and return.
     if samplegeodict is None:
         data, src = _read_pixels(src, None)
-        gd = _get_geodict_from_src(src)
+        gd = geodict_from_src(src)
         grid = Grid2D(data, gd)
         src.close()
         return grid
@@ -373,7 +365,10 @@ def read(filename, samplegeodict=None, resample=False,
         grid = Grid2D(data, gd)
 
     if resample:
-        grid = grid.interpolateToGrid(samplegeodict, method=method)
+        if interp_approach == 'scipy':
+            grid = grid.interpolateToGrid(samplegeodict, method=method)
+        else:
+            grid = grid.interpolate2(samplegeodict, method=method)
 
     src.close()
     return grid

mapio/shake.py

@@ -199,7 +199,7 @@ def readShakeFile(fileobj, adjust='bounds'):
     columns.insert(0, 'lat')
     columns.insert(0, 'lon')
     dframe = pd.read_csv(fileobj, sep='\s+', names=columns,
-                         header=None, comment='<')
+                         header=None, comment='<', dtype=np.float32)
     for field in fields:
         layers[field] = dframe[field].values.reshape(ny, nx)