diff --git a/plio/__init__.py b/plio/__init__.py index d3f342b7d48cd4d52201ac799c685be5755e2823..3e4efb6ce20e714db61902d878bc680fb51b0b27 100755 --- a/plio/__init__.py +++ b/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 - diff --git a/plio/geofuncs/__init__.py b/plio/geofuncs/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/plio/geofuncs/geofuncs.py b/plio/geofuncs/geofuncs.py new file mode 100644 index 0000000000000000000000000000000000000000..1437a58eb587fd304e365a258a95b42b98939fed --- /dev/null +++ b/plio/geofuncs/geofuncs.py @@ -0,0 +1,38 @@ +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 diff --git a/plio/io/io_gdal.py b/plio/io/io_gdal.py index b794d5d292fee90e9b309f5b3f25874e331e0d51..4afb41d8ce336df63adf9572bef1b83b7fcd7126 100644 --- a/plio/io/io_gdal.py +++ b/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)) diff --git a/plio/io/io_pvl.py b/plio/io/io_pvl.py new file mode 100644 index 0000000000000000000000000000000000000000..79c4574be08cd8554a98bf4b080b4275d50ce709 --- /dev/null +++ b/plio/io/io_pvl.py @@ -0,0 +1,20 @@ +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 diff --git a/plio/utils/utils.py b/plio/utils/utils.py index 1b97f2407fdb156262c44b39d2454890d793d4ea..761f70d28416f338b896d10f1daff6c4d7eb3ad2 100644 --- a/plio/utils/utils.py +++ b/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)