Merge pull request #11 from jlaura/master

MBR Computation

.travis.yml

@@ -45,12 +45,11 @@ install:
   - conda config --add channels jlaura
   - conda install python=$PYTHON_VERSION
   - conda install -c conda-forge gdal h5py
-  - conda install pandas sqlalchemy pyyaml networkx
+  - conda install pandas sqlalchemy pyyaml networkx affine
   - conda install -c jlaura pvl protobuf
 
   # Development installation
-  - conda install nose coverage sh anaconda-client
-  - pip install coveralls
+  - conda install pytest  pytest-cov sh anaconda-client
 
   # Straight from the menpo team
   - if [["$TRAVIS_OS_NAME" == "osx"]]; then
@@ -62,7 +61,7 @@ install:
 
 
 script:
-  - nosetests --with-coverage --cover-package=plio
+  - pytest --cov=plio 
   # After test success, package and upload to Anaconda
   - ~/miniconda/bin/python condaci.py build ./conda
 

appveyor.yml

@@ -40,7 +40,7 @@ install:
     - cmd: set PYTHONUNBUFFERED=1
 
     - cmd: conda config --set show_channel_urls true
-    - cmd: conda install --yes python=3.5 
+    - cmd: conda install --yes python=3.5
     - cmd: conda install -c pelson/channel/development --yes --quiet obvious-ci
     - cmd: conda config --add channels conda-forge
     - cmd: conda info
@@ -54,18 +54,18 @@ install:
     - cmd: conda config --add channels conda-forge
     - cmd: conda config --add channels jlaura
     - cmd: conda install --yes -c conda-forge gdal h5py
-    - cmd: conda install --yes pandas sqlalchemy pyyaml networkx
+    - cmd: conda install --yes pandas sqlalchemy pyyaml networkx affine
     - cmd: conda install --yes -c jlaura protobuf pvl
 
     # Development installation
-    - cmd: conda install --yes nose coverage
+    - cmd: conda install --yes pytest pytest-cov
     - cmd: pip install coveralls
 
 # Skip .NET project specific build phase.
 build: off
 
 test_script:
-    - cmd: nosetests --with-coverage --cover-package=plio
+    - cmd: pytest --cov=plio --ignore=plio/examples
     - "%CMD_IN_ENV% conda build conda --quiet"
 
 deploy_script:

conda/meta.yaml

@@ -23,6 +23,8 @@ requirements:
     - pandas
     - sqlalchemy
     - pyyaml
+    - affine
+    - networkx
   run:
     - python
     - setuptools
@@ -35,6 +37,8 @@ requirements:
     - pandas
     - sqlalchemy
     - pyyaml
+    - affine
+    - networkx
 
 test:
   imports:

plio/__init__.py

@@ -5,5 +5,5 @@ from . import io
 from . import date
 from . import data
 from . import examples
+from . import geofuncs
 from . import utils
-

plio/date/__init__.py

+from . import julian2ls
+from . import julian2season

plio/geofuncs/geofuncs.py

+import json
+import numpy as np
+
+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)
+    """
+    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
+
+    if minx < 0:
+        minx = 0
+    if miny < 0:
+        miny = 0
+
+    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], bufferd)
+    b_slice = intersection_to_pixels(geodata_b.latlon_to_pixel,
+                                    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 ul, ur, lr, ll
+
+def is_clockwise(vertices):
+    """
+    Returns whether a list of points describing a polygon are clockwise or counterclockwise.
+    is_clockwise(Point list) -> bool
+    Parameters
+    ----------
+    vertices : a list of points that form a single ring
+    Examples
+    --------
+    >>> is_clockwise([Point((0, 0)), Point((10, 0)), Point((0, 10))])
+    False
+    >>> is_clockwise([Point((0, 0)), Point((0, 10)), Point((10, 0))])
+    True
+    >>> v = [(-106.57798, 35.174143999999998), (-106.583412, 35.174141999999996), (-106.58417999999999, 35.174143000000001), (-106.58377999999999, 35.175542999999998), (-106.58287999999999, 35.180543), (-106.58263099999999, 35.181455), (-106.58257999999999, 35.181643000000001), (-106.58198299999999, 35.184615000000001), (-106.58148, 35.187242999999995), (-106.58127999999999, 35.188243), (-106.58138, 35.188243), (-106.58108, 35.189442999999997), (-106.58104, 35.189644000000001), (-106.58028, 35.193442999999995), (-106.580029, 35.194541000000001), (-106.57974399999999, 35.195785999999998), (-106.579475, 35.196961999999999), (-106.57922699999999, 35.198042999999998), (-106.578397, 35.201665999999996), (-106.57827999999999, 35.201642999999997), (-106.57737999999999, 35.201642999999997), (-106.57697999999999, 35.201543000000001), (-106.56436599999999, 35.200311999999997), (-106.56058, 35.199942999999998), (-106.56048, 35.197342999999996), (-106.56048, 35.195842999999996), (-106.56048, 35.194342999999996), (-106.56048, 35.193142999999999), (-106.56048, 35.191873999999999), (-106.56048, 35.191742999999995), (-106.56048, 35.190242999999995), (-106.56037999999999, 35.188642999999999), (-106.56037999999999, 35.187242999999995), (-106.56037999999999, 35.186842999999996), (-106.56037999999999, 35.186552999999996), (-106.56037999999999, 35.185842999999998), (-106.56037999999999, 35.184443000000002), (-106.56037999999999, 35.182943000000002), (-106.56037999999999, 35.181342999999998), (-106.56037999999999, 35.180433000000001), (-106.56037999999999, 35.179943000000002), (-106.56037999999999, 35.178542999999998), (-106.56037999999999, 35.177790999999999), (-106.56037999999999, 35.177143999999998), (-106.56037999999999, 35.175643999999998), (-106.56037999999999, 35.174444000000001), (-106.56037999999999, 35.174043999999995), (-106.560526, 35.174043999999995), (-106.56478, 35.174043999999995), (-106.56627999999999, 35.174143999999998), (-106.566541, 35.174144999999996), (-106.569023, 35.174157000000001), (-106.56917199999999, 35.174157999999998), (-106.56938, 35.174143999999998), (-106.57061499999999, 35.174143999999998), (-106.57097999999999, 35.174143999999998), (-106.57679999999999, 35.174143999999998), (-106.57798, 35.174143999999998)]
+    >>> is_clockwise(v)
+    True
+    """
+    if len(vertices) < 3:
+        return True
+    area = 0.0
+    ax, ay = vertices[0]
+    for bx, by in vertices[1:]:
+        area += ax * by - ay * bx
+        ax, ay = bx, by
+    bx, by = vertices[0]
+    area += ax * by - ay * bx
+    return area < 0.0

plio/geofuncs/tests/test_geofuncs.py

+import os
+import sys
+sys.path.insert(0, os.path.abspath('..'))
+
+import pytest
+
+from plio.geofuncs import geofuncs
+
+@pytest.fixture
+def nevada():
+    return [[-114.04392,40.68928],[-114.04558,40.4958],[-114.04619,40.30302],[-114.04644,40.09896],[-114.04658,39.99994],[-114.04727,39.75817],[-114.04757,39.61018],[-114.0473,39.45715],[-114.04779,39.36296],[-114.04841,39.23851],[-114.04885,39.08777],[-114.04833,38.90545],[-114.04916,38.75165],[-114.04992,38.55049],[-114.04997,38.20495],[-114.05013,37.95499],[-114.04939,37.77873],[-114.05198,37.70735],[-114.05264,37.47222],[-114.05187,37.13439],[-114.0506,37.0004],[-114.0506,36.99997],[-114.05014,36.817],[-114.04736,36.60322],[-114.04338,36.37619],[-114.04404,36.21464],[-114.1139,36.09833],[-114.22646,36.01461],[-114.32346,36.10119],[-114.51122,36.15058],[-114.6729,36.11546],[-114.73513,36.05493],[-114.74365,35.98542],[-114.70883,35.9167],[-114.67489,35.86436],[-114.70415,35.81412],[-114.69704,35.73579],[-114.68702,35.66942],[-114.65449,35.60517],[-114.66076,35.5417],[-114.6768,35.49125],[-114.61121,35.37012],[-114.58031,35.21811],[-114.57354,35.14231],[-114.63064,35.11791],[-114.60899,35.07971],[-114.63423,35.00332],[-114.63349,35.00186],[-114.63361,35.00195],[-114.82052,35.15341],[-115.11622,35.38796],[-115.36992,35.59033],[-115.65233,35.81231],[-115.89512,36.0018],[-116.08072,36.14577],[-116.37528,36.37205],[-116.87227,36.75057],[-117.31883,37.08441],[-117.79563,37.43715],[-118.04392,37.6185],[-118.22972,37.75309],[-118.51722,37.96065],[-119.00097,38.30368],[-119.43506,38.60904],[-119.76041,38.83427],[-119.9748,38.98156],[-120.00608,39.37557],[-120.0015,39.57782],[-120.00049,39.79567],[-119.99733,40.08934],[-119.99567,40.39719],[-119.99926,40.86934],[-120.00002,41.26742],[-119.99919,41.97905],[-119.99917,41.99454],[-119.99917,41.99484],[-119.90622,41.9972],[-119.80128,41.99746],[-119.70479,41.99621],[-119.61469,41.99594],[-119.48157,41.99572],[-119.36302,41.99428],[-119.25103,41.99384],[-119.0022,41.99375],[-118.7824,41.99259],[-118.54194,41.99478],[-118.42657,41.99586],[-118.3477,41.99629],[-118.19842,41.99699],[-118.09918,41.9975],[-118.06386,41.99767],[-117.99278,41.99804],[-117.90359,41.99817],[-117.79003,41.99842],[-117.71943,41.99819],[-117.62441,41.99836],[-117.51839,41.99915],[-117.46483,41.9994],[-117.33055,41.9996],[-117.18204,42.00038],[-117.02916,42.00018],[-117.02622,42.00025],[-117.02234,41.99989],[-116.861,41.99876],[-116.70972,41.99812],[-116.55307,41.99766],[-116.48554,41.99686],[-116.37594,41.99632],[-116.20623,41.99768],[-116.12593,41.99765],[-116.01199,41.99804],[-115.98839,41.99855],[-115.84911,41.99677],[-115.69522,41.99701],[-115.52677,41.99676],[-115.3968,41.99648],[-115.31388,41.9961],[-115.25015,41.99615],[-115.12254,41.99607],[-114.91817,41.99977],[-114.76097,41.99991],[-114.6492,41.9962],[-114.61271,41.99503],[-114.49854,41.9946],[-114.39727,41.99501],[-114.30408,41.99437],[-114.18492,41.99407],[-114.04537,41.99372],[-114.04172,41.99372],[-114.0398,41.89425],[-114.04055,41.59062],[-114.04061,41.36],[-114.04195,41.05548],[-114.04375,40.76026],[-114.04391,40.68985]]
+
+def test_find_four_corners(nevada):
+    corners = geofuncs.find_four_corners(nevada)
+    assert corners[0] == [-119.99917, 41.99484]
+    assert corners[3] == [-114.63349, 35.00186]
+    assert len(corners) == 4
+
+def test_find_corners(nevada):
+    corners = geofuncs.find_corners(nevada)
+    assert len(corners) == 8
+
+    corners = geofuncs.find_corners(nevada, threshold=100)
+    assert len(corners) == 5

plio/io/__init__.py

+from . import io_autocnetgraph
+from . import io_controlnetwork
+from . import io_db
+from . import io_gdal
+from . import io_hdf
+from . import io_json
+from . import io_krc
+from . import io_pvl
+from . import io_spectral_profiler
+from . import io_tes
+from . import io_yaml
+from . import isis_serial_number

plio/io/io_gdal.py

+import json
+import math
 import os
 import warnings
 
+import affine
 import gdal
 import numpy as np
 import osr
@@ -8,6 +11,7 @@ import pvl
 from osgeo import ogr
 
 from plio.io import extract_metadata
+from plio.geofuncs import geofuncs
 from plio.utils.utils import find_in_dict
 
 gdal.UseExceptions()
@@ -27,7 +31,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()):
@@ -173,22 +178,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):
@@ -201,6 +215,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 geofuncs.is_clockwise(json.loads(self.footprint.ExportToJson())['coordinates'][0][0])
+        else:
+            return True
+
     @property
     def spatial_reference(self):
         if not getattr(self, '_srs', None):
@@ -396,16 +417,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):
@@ -424,11 +437,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 * (lon, lat))
+        return px, py
 
     def read_array(self, band=1, pixels=None, dtype='float32'):
         """
@@ -458,25 +469,34 @@ 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 compute_overlap(self, geodata, **kwargs):
+        return geofuncs.compute_overlap(self, geodata, **kwargs)
+
 
 def array_to_raster(array, file_name, projection=None,
                     geotransform=None, outformat='GTiff',
@@ -549,7 +569,7 @@ def array_to_raster(array, file_name, projection=None,
 
 
 def match_rasters(match_to, match_from, destination,
-                  resampling_method='GRA_Bilinear'):
+                  resampling_method='GRA_Bilinear', ndv=0):
     """
     Match a source raster to a match raster, including resolution and extent.
 
@@ -590,7 +610,9 @@ def match_rasters(match_to, match_from, destination,
 
     dst = gdal.GetDriverByName('GTiff').Create(destination, width, height, 1,
                                                gdalconst.GDT_Float64)
+
     dst.SetGeoTransform(match_to_gt)
     dst.SetProjection(match_to_srs)
+    dst.GetRasterBand(1).SetNoDataValue(ndv)
 
     gdal.ReprojectImage(match_from.dataset, dst, None, None, getattr(gdalconst, resampling_method))

plio/io/io_pvl.py

+import pvl
+from plio.utils.utils import find_in_dict
+
+def extract_keywords(header, *args):
+    """
+    For a given header, find all of the keys and return an unnested dict.
+    """
+    try:
+        header = pvl.load(header)
+    except:
+        header = pvl.loads(header)
+
+    res = {}
+    # Iterate through all of the requested keys
+    for a in args:
+        try:
+            res[a] = find_in_dict(a)
+        except:
+            res[a] = None
+    return res

plio/io/tests/test_io_gdal.py

@@ -189,7 +189,7 @@ class TestWriter(unittest.TestCase):
         no_data_value = 0.0
         io_gdal.array_to_raster(self.arr, 'test.tif', ndv=no_data_value)
         dataset = io_gdal.GeoDataset('test.tif')
-        self.assertEqual(dataset.no_data_value, no_data_value) 
+        self.assertEqual(dataset.no_data_value, no_data_value)
 
     def test_with_projection(self):
         wktsrs = """PROJCS["Moon2000_Mercator180",
@@ -226,6 +226,3 @@ class TestWriter(unittest.TestCase):
             os.remove('test.tif')
         except:
             pass
-
-
-

plio/utils/utils.py

@@ -18,6 +18,16 @@ def create_dir(basedir=''):
     """
     return tempfile.mkdtemp(dir=basedir)
 
+def check_file_exists(fpath):
+    #Ensure that the file exists at the PATH specified
+    if os.path.isfile(fpath) == False:
+        error_msg = "Unable to find file: {}\n".format(fpath)
+        try:
+            logging.error(error_msg)
+        except:
+            print(error_msg)
+        return
+    return True
 
 def delete_dir(dir):
     """
@@ -153,4 +163,4 @@ def xstr(s):
     """
     if s is None:
         return ''
-    return str(s)
\ No newline at end of file
+    return str(s)