# -*- coding: utf-8 -*- # py_tools_ds - A collection of geospatial data analysis tools that simplify standard # operations when handling geospatial raster and vector data as well as projections. # # Copyright (C) 2016-2021 # - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de) # - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam, # Germany (https://www.gfz-potsdam.de/) # # This software was developed within the context of the GeoMultiSens project funded # by the German Federal Ministry of Education and Research # (project grant code: 01 IS 14 010 A-C). # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. # You may obtain a copy of the License at # # http://www.apache.org/licenses/LICENSE-2.0 # # Unless required by applicable law or agreed to in writing, software # distributed under the License is distributed on an "AS IS" BASIS, # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. import math import warnings import numpy as np from typing import Union # noqa F401 # flake8 issue from geopandas import GeoDataFrame from shapely.geometry import shape, Polygon, box from shapely.geometry import MultiPolygon # noqa F401 # flake8 issue from ..coord_trafo import mapYX2imYX from ..coord_grid import find_nearest_grid_coord from ..coord_calc import get_corner_coordinates __author__ = "Daniel Scheffler" def get_overlap_polygon(poly1, poly2, v=False): """Return a dict with the overlap of two shapely.Polygon() objects, the overlap percentage and the overlap area. :param poly1: first shapely.Polygon() object :param poly2: second shapely.Polygon() object :param v: verbose mode :return: overlap polygon as shapely.Polygon() object :return: overlap percentage as float value [%] :return: area of overlap polygon """ # compute overlap polygon overlap_poly = poly1.intersection(poly2) if overlap_poly.geom_type == 'GeometryCollection': overlap_poly = overlap_poly.buffer(0) # converts 'GeometryCollection' to 'MultiPolygon' if not overlap_poly.is_empty: # check if output is MultiPolygon or GeometryCollection -> if yes, convert to Polygon if overlap_poly.geom_type == 'MultiPolygon': overlap_poly = fill_holes_within_poly(overlap_poly) assert overlap_poly.geom_type == 'Polygon', \ "get_overlap_polygon() did not return geometry type 'Polygon' but %s." % overlap_poly.geom_type overlap_percentage = 100 * shape(overlap_poly).area / shape(poly2).area if v: print('%.2f percent of the image to be shifted is covered by the reference image.' % overlap_percentage) # pragma: no cover return {'overlap poly': overlap_poly, 'overlap percentage': overlap_percentage, 'overlap area': overlap_poly.area} else: return {'overlap poly': None, 'overlap percentage': 0, 'overlap area': 0} def get_footprint_polygon(CornerLonLat, fix_invalid=False): """Convert a list of coordinates into a shapely polygon object. :param CornerLonLat: a list of coordinate tuples like [[lon,lat], [lon. lat], ..] in clockwise or counter clockwise order :param fix_invalid: fix invalid output polygon by returning its convex hull (sometimes this can be different) :return: a shapely.Polygon() object """ if fix_invalid: with warnings.catch_warnings(): warnings.simplefilter('ignore') # FIXME not working outpoly = Polygon(CornerLonLat) if not outpoly.is_valid: outpoly = outpoly.convex_hull else: outpoly = Polygon(CornerLonLat) assert outpoly.is_valid, 'The given coordinates result in an invalid polygon. Check coordinate order.' \ 'Got coordinates %s.' % CornerLonLat return outpoly def get_bounds_polygon(gt, rows, cols): """Get a polygon representing the outer bounds of an image. :param gt: GDAL geotransform :param rows: number of rows :param cols: number of columns :return: """ return get_footprint_polygon(get_corner_coordinates(gt=gt, cols=cols, rows=rows)) def get_smallest_boxImYX_that_contains_boxMapYX(box_mapYX, gt_im, tolerance_ndigits=5): """Return image coordinates of the smallest box at the given coordinate grid that contains the given map coords box. :param box_mapYX: input box coordinates as YX-tuples :param gt_im: geotransform of input box :param tolerance_ndigits: tolerance to avoid that output image coordinates are rounded to next integer although they have been very close to an integer before (this avoids float rounding issues) -> tolerance is given as number of decimal digits of an image coordinate :return: """ xmin, ymin, xmax, ymax = Polygon([(i[1], i[0]) for i in box_mapYX]).bounds # map-bounds box_mapYX (ymin, xmin), (ymax, xmax) = mapYX2imYX([ymin, xmin], gt_im), mapYX2imYX([ymax, xmax], gt_im) # image coord bounds # round min coords off and max coords on but tolerate differences below n decimal digits as the integer itself xmin, ymin, xmax, ymax = np.round([xmin, ymin, xmax, ymax], tolerance_ndigits) xmin, ymin, xmax, ymax = math.floor(xmin), math.ceil(ymin), math.ceil(xmax), math.floor(ymax) return (ymax, xmin), (ymax, xmax), (ymin, xmax), (ymin, xmin) # UL_YX,UR_YX,LR_YX,LL_YX def get_largest_onGridPoly_within_poly(outerPoly, gt, rows, cols): oP_xmin, oP_ymin, oP_xmax, oP_ymax = outerPoly.bounds xmin, ymax = find_nearest_grid_coord((oP_xmin, oP_ymax), gt, rows, cols, direction='SE') xmax, ymin = find_nearest_grid_coord((oP_xmax, oP_ymin), gt, rows, cols, direction='NW') return box(xmin, ymin, xmax, ymax) def get_smallest_shapelyImPolyOnGrid_that_contains_shapelyImPoly(shapelyPoly): """Return the smallest box that matches the coordinate grid of the given geotransform. The returned shapely polygon contains image coordinates.""" xmin, ymin, xmax, ymax = shapelyPoly.bounds # image_coords-bounds # round min coords off and max coords on out_poly = box(math.floor(xmin), math.floor(ymin), math.ceil(xmax), math.ceil(ymax)) return out_poly def find_line_intersection_point(line1, line2): xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0]) ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1]) def det(a, b): return a[0] * b[1] - a[1] * b[0] div = det(xdiff, ydiff) if div == 0: return None, None d = (det(*line1), det(*line2)) x = det(d, xdiff) / div y = det(d, ydiff) / div return x, y def polyVertices_outside_poly(inner_poly, outer_poly, tolerance=0): """Check if a shapely polygon (inner_poly) contains vertices that are outside of another polygon (outer_poly). :param inner_poly: the polygon with the vertices to check :param outer_poly: the polygon where all vertices have to be inside :param tolerance: tolerance of the decision """ return not outer_poly.buffer(tolerance).contains(inner_poly) def fill_holes_within_poly(poly): # type: (Union[Polygon, MultiPolygon]) -> Polygon """Fill the holes within a shapely Polygon or MultiPolygon and return a Polygon with only the outer boundary. :param poly: , shapely.geometry.GeometryCollection> :return: """ if poly.geom_type not in ['Polygon', 'MultiPolygon']: raise ValueError("Unexpected geometry type %s." % poly.geom_type) if poly.geom_type == 'Polygon': # return only the exterior polygon filled_poly = Polygon(poly.exterior) else: # 'MultiPolygon' gdf = GeoDataFrame(columns=['geometry']) gdf['geometry'] = poly.geoms # get the area of each polygon of the multipolygon EXCLUDING the gaps in it # noinspection PyTypeChecker gdf['area_filled'] = gdf.apply( lambda GDF_row: Polygon(np.swapaxes(np.array(GDF_row.geometry.exterior.coords.xy), 0, 1)).area, axis=1) largest_poly_filled = gdf.loc[gdf['area_filled'].astype(np.float64).idxmax()]['geometry'] # noinspection PyTypeChecker gdf['within_or_equal'] = gdf.apply( lambda GDF_row: GDF_row.geometry.within(largest_poly_filled.buffer(1e-5)) or GDF_row.geometry.equals(largest_poly_filled), axis=1) if False in gdf.within_or_equal.values: n_disjunct_polys = int(np.sum(~gdf.within_or_equal)) warnings.warn(RuntimeWarning('The given MultiPolygon contains %d disjunct polygone(s) outside of the ' 'largest polygone. fill_holes_within_poly() will only return the largest ' 'polygone as a filled version.' % n_disjunct_polys)) # return the outer boundary of the largest polygon filled_poly = Polygon(np.swapaxes(np.array(largest_poly_filled.exterior.coords.xy), 0, 1)) return filled_poly