Addition for affine transformations

plio/geofuncs/geofuncs.py

+import json
 import numpy as np
 
-def intersection_to_pixels(latlon_to_pixel, ul, lr, xmax, ymax, buffer=300):
+def intersection_to_pixels(inverse_affine, ul, ur, lr, ll):
     """
+    Given an inverse affine transformation, take four bounding coordinates
+    in lat/lon space and convert them to a minimum bounding rectangle in pixel
+    space.
 
+    inverse_affine : object
+                     An Affine transformation object
+
+    ul : tuple
+         Upper Left coordinate in the form (x,y)
+
+    ur : tuple
+         Upper Right coordinate in the form (x,y)
+
+    lr : tuple
+         Lower Right coordinate in the form (x,y)
+
+    ll : tuple
+         Lower Left coordinate in the form (x,y)
     """
-    x_start, y_start = latlon_to_pixel(ul[1], ul[0])
-    x_stop, y_stop = latlon_to_pixel(lr[1], lr[0])
-
-    if x_start > x_stop:
-        x_start, x_stop = x_stop, x_start
-    if y_start > y_stop:
-        y_start, y_stop = y_stop, y_start
-
-    if x_start < buffer:
-        x_start = 0
-    if xmax - x_stop < buffer:
-        x_stop = xmax
-
-    if y_start < buffer:
-        y_start = 0
-    if ymax - y_stop < buffer:
-        y_stop = ymax
-    return np.s_[y_start:y_stop, x_start:x_stop]
-
-def estimate_mbr(geodata_a, geodata_b):
+    minx = np.inf
+    maxx = -np.inf
+    miny = np.inf
+    maxy = -np.inf
+
+    for c in [ul, ur, lr, ll]:
+        px, py = map(int, inverse_affine * (c[0], c[1]))
+
+        if px < minx:
+            minx = px
+        if px > maxx:
+            maxx = px
+
+        if py < miny:
+            miny = py
+        if py > maxy:
+            maxy = py
+
+    return minx, maxx, miny, maxy
+
+def compute_overlap(geodata_a, geodata_b):
+    p1 = geodata_a.footprint
+    p2 = geodata_b.footprint
+    intersection = json.loads(p1.Intersection(p2).ExportToJson())['coordinates'][0]
+
+    ul, ur, lr, ll = find_four_corners(intersection)
+
+    a_intersection = intersection_to_pixels(geodata_a.inverse_affine, ul, ur, lr, ll)
+    b_intersection = intersection_to_pixels(geodata_b.inverse_affine, ul, ur, lr, ll)
+
+    return (ul, ur, lr, ll), a_intersection, b_intersection
+
+def estimate_mbr(geodata_a, geodata_b, bufferd=50):
     p1 = geodata_a.footprint
     p2 = geodata_b.footprint
     minx, maxx, miny, maxy = p1.Intersection(p2).GetEnvelope()
     ul = (minx, maxy)
     lr = (maxx, miny)
-
     a_slice = intersection_to_pixels(geodata_a.latlon_to_pixel,
-                                    ul, lr, *geodata_a.xy_extent[1])
+                                    ul, lr, *geodata_a.xy_extent[1], bufferd)
     b_slice = intersection_to_pixels(geodata_b.latlon_to_pixel,
-                                    ul, lr, *geodata_b.xy_extent[1])
+                                    ul, lr, *geodata_b.xy_extent[1], bufferd)
+
+    return a_slice, b_slice
+
+
+def find_corners(coords, threshold=120):
+    """
+    Given a list of coordinates in the form [(x, y), (x1, y1), ..., (xn, yn)],
+    attempt to find corners by identifying all angles < the threshold.  For a
+    line segment composed of 3 points, the angle between ab and ac is 180 if the
+    segments are colinear.  If they are less than threshold, the line segments
+    are considered to be a corner.
+
+    Parameters
+    ----------
+    coords : list
+             of coordinates
+
+    threshold : numeric
+                Angle under which a corner is identified, Default: 120
+    """
+    corners = [coords[0]]
+    for i, a in enumerate(coords[:-2]):
+        a = np.asarray(a)
+        b = np.asarray(coords[i+1])
+        c = np.asarray(coords[i+2])
+        ba = a - b
+        bc = c - b
+        angle = np.arccos(np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc)))
+        if np.degrees(angle) < threshold:
+            corners.append(b.tolist())
+
+    return corners
+
+def find_four_corners(coords, threshold=120):
+    """
+    Find four corners in a polygon by making a call to the find_corners
+    function and using the first four corners.
+
+    Parameters
+    ----------
+    coords: list
+            of coordinates
+
+    threshold : numeric
+                Anfle under which a corner is identified, Default: 120
+
+    See Also
+    --------
+    plio.geofuncs.geofuncs.find_corners
+    """
+    corners = find_corners(coords, threshold)
+
+    corners.sort(key = lambda x:x[1])
+    upper = corners[2:]
+    lower = corners[:2]
+
+
+    key = lambda x:x[0]
+    ul = min(upper, key=key)
+    ur = max(upper, key=key)
+    lr = max(lower, key=key)
+    ll = min(lower, key=key)
 
-    return a_slice, b_slice                                
+    return ul, ur, lr, ll

plio/io/io_gdal.py

+import json
+import math
 import os
 import warnings
 
+import affine
 import gdal
 import numpy as np
 import osr
@@ -11,6 +14,8 @@ from plio.io import extract_metadata
 from plio.geofuncs import geofuncs
 from plio.utils.utils import find_in_dict
 
+from libpysal.cg import is_clockwise
+
 gdal.UseExceptions()
 
 NP2GDAL_CONVERSION = {
@@ -28,7 +33,8 @@ NP2GDAL_CONVERSION = {
 
 GDAL2NP_CONVERSION = {}
 
-DEFAULT_PROJECTIONS = {'mars':'GEOGCS["Mars 2000",DATUM["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,169.89444722361179]],PRIMEM["Greenwich",0],UNIT["Decimal_Degree",0.0174532925199433]]'}
+DEFAULT_PROJECTIONS = {'mars':'GEOGCS["Mars 2000",DATUM["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,169.89444722361179]],PRIMEM["Greenwich",0],UNIT["Decimal_Degree",0.0174532925199433]]',
+                       'moon':'GEOGCS["Moon 2000",DATUM["D_Moon_2000",SPHEROID["Moon_2000_IAU_IAG",1737400.0,0.0]],PRIMEM["Greenwich",0],UNIT["Decimal_Degree",0.0174532925199433]]'}
 DEFAULT_RADII = {'mars': 3396190.0}
 
 for k, v in iter(NP2GDAL_CONVERSION.items()):
@@ -174,22 +180,31 @@ class GeoDataset(object):
         """
         if not getattr(self, '_geotransform', None):
             if self.footprint:
-                # This is an ISIS3 cube that does not report a valid geotransform
-                ul, ll, lr, ur = self.latlon_extent
-                xs, ys = self.raster_size
-                xres = abs(ul[0] - ur[0]) / xs
-                yres = abs(ul[1] - ll[1]) / ys
-                self._geotransform = [ul[0], xres, 0, ul[1], 0, -yres]
-
-                if ul[1] > ll[1]:
-                    # Image is South-North instead of North-South
-                    self._geotransform[-1] *= -1  # Lat increases with image length
-                    self._geotransform[3] = ll[1] # Origin is at the LL
-
+                coords = json.loads(self.footprint.ExportToJson())['coordinates'][0][0]
+                ul, ur, lr, ll = geofuncs.find_four_corners(coords)
+
+                xsize, ysize = self.raster_size
+                xstart = ul[0]
+                xscale = (ur[0] - xstart) / xsize
+                xskew = (ll[0] - xstart) / ysize
+                ystart = ul[1]
+                yskew = (ur[1] - ystart) / xsize
+                yscale = (ll[1] - ystart) / ysize
+                self._geotransform = [xstart, xscale, xskew, ystart, yskew, yscale]
             else:
                 self._geotransform = self.dataset.GetGeoTransform()
         return self._geotransform
 
+    @property
+    def forward_affine(self):
+        self._fa = affine.Affine.from_gdal(*self.geotransform)
+        return self._fa
+
+    @property
+    def inverse_affine(self):
+        self._ia = ~self.forward_affine
+        return self._ia
+
     @property
     def standard_parallels(self):
         if not getattr(self, '_standard_parallels', None):
@@ -202,6 +217,13 @@ class GeoDataset(object):
             self._unit_type = self.dataset.GetRasterBand(1).GetUnitType()
         return self._unit_type
 
+    @property
+    def north_up(self):
+        if self.footprint:
+            return is_clockwise(json.loads(self.footprint.ExportToJson())['coordinates'][0][0])
+        else:
+            return True
+
     @property
     def spatial_reference(self):
         if not getattr(self, '_srs', None):
@@ -397,16 +419,8 @@ class GeoDataset(object):
                    (Latitude, Longitude) corresponding to the given (x,y).
 
         """
-        try:
-            geotransform = self.geotransform
-            x = geotransform[0] + (x * geotransform[1]) + (y * geotransform[2])
-            y = geotransform[3] + (x * geotransform[4]) + (y * geotransform[5])
-            lon, lat, _ = self.coordinate_transformation.TransformPoint(x, y)
-        except:
-            lat = lon = None
-            warnings.warn('Unable to compute pixel to geographic conversion without '
-                          'projection information for {}'.format(self.base_name))
-
+        lon, lat = self.forward_affine * (x,y)
+        lon, lat, _ = self.coordinate_transformation.TransformPoint(lon, lat)
         return lat, lon
 
     def latlon_to_pixel(self, lat, lon):
@@ -425,11 +439,9 @@ class GeoDataset(object):
                (Sample, line) position corresponding to the given (latitude, longitude).
 
         """
-        geotransform = self.geotransform
-        upperlat, upperlon, _ = self.inverse_coordinate_transformation.TransformPoint(lon, lat)
-        x = (upperlat - geotransform[0]) / geotransform[1]
-        y = (upperlon - geotransform[3]) / geotransform[5]
-        return int(x), int(y)
+        lon, lat, _ = self.inverse_coordinate_transformation.TransformPoint(lon, lat)
+        px, py = map(int, self.inverse_affine * (lat, lon))
+        return px, py
 
     def read_array(self, band=1, pixels=None, dtype='float32'):
         """
@@ -459,27 +471,33 @@ class GeoDataset(object):
 
         if not pixels:
             array = band.ReadAsArray().astype(dtype)
+            if self.north_up == False:
+                array = np.flipud(array)
         else:
             # Check that the read start is not outside of the image
-            xstart, ystart, xextent, yextent = pixels
+            xstart, ystart, xcount, ycount = pixels
+            xmax, ymax = map(int, self.xy_extent[1])
+
+            # If the image is south up, flip the roi
+            if self.north_up == False:
+                ystart = ymax - ystart - ycount
             if xstart < 0:
                 xstart = 0
 
             if ystart < 0:
                 ystart = 0
 
-            xmax, ymax = map(int, self.xy_extent[1])
-            if xstart + pixels[2] > xmax:
-                xextent = xmax - xstart
+            if xstart + xcount > xmax:
+                xcount = xmax - xstart
 
-            if ystart + pixels[3] > ymax:
-                yextent = ymax - ystart
+            if ystart +  ycount > ymax:
+                ycount = ymax - ystart
 
-            array = band.ReadAsArray(xstart, ystart, xextent, yextent).astype(dtype)
+            array = band.ReadAsArray(xstart, ystart, xcount, ycount).astype(dtype)
         return array
 
-    def estimate_mbr(self, geodata):
-        return geofuncs.estimate_mbr(self, geodata)
+    def compute_overlap(self, geodata, **kwargs):
+        return geofuncs.compute_overlap(self, geodata, **kwargs)
 
 
 def array_to_raster(array, file_name, projection=None,