Commit 28948c0f authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Fixed warning within output validation of warp_ndarray in case out_bounds_prj...


Fixed warning within output validation of  warp_ndarray in case out_bounds_prj is provided in a different projection.
Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent 2d893275
......@@ -2,6 +2,12 @@
History
=======
0.17.8 (2021-08-11)
-------------------
* Fixed an issue where the output polygon of raster2polygon could contain vertices outside of the input array.
0.17.7 (2021-08-09)
-------------------
......
......@@ -29,6 +29,7 @@ from osgeo import gdal, osr, ogr
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...processing.progress_mon import ProgressBar, Timer
from ..raster.reproject import warp_ndarray
from ..vector.topology import get_bounds_polygon, polyVertices_outside_poly, get_overlap_polygon
def raster2polygon(array, gt, prj, DN2extract=1, exact=True, maxfeatCount=None,
......@@ -36,10 +37,10 @@ def raster2polygon(array, gt, prj, DN2extract=1, exact=True, maxfeatCount=None,
"""Calculate a footprint polygon for the given array.
:param array: 2D numpy array
:param gt:
:param prj:
:param gt: GDAL GeoTransform
:param prj: projection as WKT string, 'EPSG:1234' or <EPSG_int>
:param DN2extract: <int, float> pixel value to create polygons for
:param exact:
:param exact: whether to compute the exact footprint polygon or a simplified one for speed
:param maxfeatCount: <int> the maximum expected number of polygons. If more polygons are found, every further
processing is cancelled and a RunTimeError is raised.
:param timeout: breaks the process after a given time in seconds
......@@ -48,8 +49,10 @@ def raster2polygon(array, gt, prj, DN2extract=1, exact=True, maxfeatCount=None,
:return:
"""
assert array.ndim == 2, "Only 2D arrays are supported. Got a %sD array." % array.ndim
gt_orig = gt
shape_orig = array.shape
# downsample input array in case is has more than 1e8 pixels to prevent crash
# downsample input array in case it has more than 1e8 pixels to prevent crash
if not exact and array.size > 1e8: # 10000 x 10000 px
# downsample with nearest neighbour
zoom_factor = (8000 * 8000 / array.size) ** 0.5
......@@ -144,4 +147,13 @@ def raster2polygon(array, gt, prj, DN2extract=1, exact=True, maxfeatCount=None,
del mem_ds, mem_layer
shplyPoly = GDF.loc[1, 'geometry']
# the downsampling in case exact=False may cause vertices of shplyPoly to be outside the input array bounds
# -> clip shplyPoly with bounds_poly in that case
if not exact:
bounds_poly = get_bounds_polygon(gt_orig, *shape_orig)
if polyVertices_outside_poly(shplyPoly, bounds_poly, 1e-5):
shplyPoly = get_overlap_polygon(shplyPoly, bounds_poly)['overlap poly']
return shplyPoly
......@@ -27,10 +27,11 @@ import numpy as np
from typing import Union # noqa F401 # flake8 issue
from geopandas import GeoDataFrame
from shapely.geometry import shape, Polygon, box, Point
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"
......@@ -90,6 +91,17 @@ def get_footprint_polygon(CornerLonLat, fix_invalid=False):
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.
......@@ -145,23 +157,13 @@ def find_line_intersection_point(line1, line2):
def polyVertices_outside_poly(inner_poly, outer_poly, tolerance=0):
"""Check if a shapely polygon (inner_poly) contains vertices that do not intersect another polygon (outer_poly).
"""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
"""
if inner_poly.within(outer_poly.buffer(tolerance)):
# all vertices are inside outer_poly
return False
elif inner_poly.intersects(outer_poly.buffer(tolerance)):
# check if all vertices intersect with outer poly
GDF = GeoDataFrame(np.swapaxes(np.array(inner_poly.exterior.coords.xy), 0, 1), columns=['X', 'Y'])
# noinspection PyTypeChecker
return False in GDF.apply(lambda GDF_row: Point(GDF_row.X, GDF_row.Y).intersects(outer_poly), axis=1).values
else:
# inner_poly does not intersect out_poly -> all vertices are outside
return True
return not outer_poly.buffer(tolerance).contains(inner_poly)
def fill_holes_within_poly(poly):
......
......@@ -19,5 +19,5 @@
# You should have received a copy of the GNU Lesser General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
__version__ = '0.17.7'
__versionalias__ = '20210809_01'
__version__ = '0.17.8'
__versionalias__ = '20210811_01'
......@@ -32,7 +32,7 @@ Tests for `py_tools_ds.io.raster.gdal` module.
from unittest import TestCase, main
import numpy as np
from gdal import Dataset
from osgeo.gdal import Dataset
from pandas import DataFrame
from py_tools_ds.io.raster.gdal import get_GDAL_ds_inmem, get_GDAL_driverList
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment