Adds the ability to estimate a minimum bounding rectangle between two geo datasets

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/geofuncs/geofuncs.py

+import numpy as np
+
+def intersection_to_pixels(latlon_to_pixel, ul, lr, xmax, ymax, buffer=300):
+    """
+
+    """
+    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):
+    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])
+    b_slice = intersection_to_pixels(geodata_b.latlon_to_pixel,
+                                    ul, lr, *geodata_b.xy_extent[1])
+
+    return a_slice, b_slice                                

plio/io/io_gdal.py

@@ -8,6 +8,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()
@@ -477,6 +478,9 @@ class GeoDataset(object):
             array = band.ReadAsArray(xstart, ystart, xextent, yextent).astype(dtype)
         return array
 
+    def estimate_mbr(self, geodata):
+        return geofuncs.estimate_mbr(self, geodata)
+
 
 def array_to_raster(array, file_name, projection=None,
                     geotransform=None, outformat='GTiff',
@@ -549,7 +553,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 +594,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/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)