diff --git a/plio/geofuncs/geofuncs.py b/plio/geofuncs/geofuncs.py index 1437a58eb587fd304e365a258a95b42b98939fed..94e7d07af4c41de76cbf78c5f09d51d5ae7362bc 100644 --- a/plio/geofuncs/geofuncs.py +++ b/plio/geofuncs/geofuncs.py @@ -1,38 +1,130 @@ +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 diff --git a/plio/io/io_gdal.py b/plio/io/io_gdal.py index 4afb41d8ce336df63adf9572bef1b83b7fcd7126..d86f8341aee72c9b935db80d42fcdec6427fa514 100644 --- a/plio/io/io_gdal.py +++ b/plio/io/io_gdal.py @@ -1,6 +1,9 @@ +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,