Commit 6df976a1 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

bugfix for raster2polygon

compatibility.python.exceptions:
- added FileNotFoundError

geo.raster.conversion:
- raster2polygon:
    - added keywords 'DN2extract' and 'maxfeatCount'
    - added dissolving of multiple little polygons within result vector layer to one multipart polygon

geo.vector.topology:
- added polyVertices_outside_poly()

io.raster.GeoArray:
- footprint_poly: added assertion and modified call of raster2polygon
- added footprint_poly setter

processing.process_mon:
- revised is_timed_out()
parent 750f15f0
......@@ -3,5 +3,9 @@ __author__ = "Daniel Scheffler"
class TimeoutError(OSError):
""" Timeout expired. """
pass
class FileNotFoundError(OSError):
""" Timeout expired. """
pass
\ No newline at end of file
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"
from six import PY3
from shapely.wkb import loads
import numpy as np
......@@ -16,16 +16,19 @@ except ImportError:
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...processing.progress_mon import printProgress, is_timed_out
from ...compatibility.python.exceptions import TimeoutError
from ...compatibility.python.exceptions import TimeoutError as TimeoutError_comp
def raster2polygon(array_or_GeoArray, gt=None, prj=None, exact=True, timeout=None, progress=True, q=False):
def raster2polygon(array_or_GeoArray, gt=None, prj=None, DN2extract=1, exact=True, maxfeatCount=None,
timeout=None, progress=True, q=False):
"""Calculates a footprint polygon for the given array or GeoArray.
:param array_or_GeoArray:
:param gt:
:param prj:
:param DN2extract:
:param exact:
:param maxfeatCount:
:param timeout: breaks the process after a given time in seconds
:param progress: show progress bars (default: True)
:param q: quiet mode (default: False)
......@@ -51,27 +54,63 @@ def raster2polygon(array_or_GeoArray, gt=None, prj=None, exact=True, timeout=Non
# set callback
progressBar = lambda percent01, message, user_data: \
printProgress(percent01 * 100, prefix='Polygonize progress ', suffix='Complete', barLength=50, timeout=3)
timeout_callback = lambda percent01, message, user_data: is_timed_out(3)
printProgress(percent01 * 100, prefix='Polygonize progress ', suffix='Complete', barLength=50, timeout=timeout)
timeout_callback = lambda percent01, message, user_data: is_timed_out(timeout)
callback = progressBar if progress and not q else timeout_callback if timeout else None
# run the algorithm
result = gdal.Polygonize(src_band, src_band.GetMaskBand(), mem_layer, 0, ["8CONNECTED=8"] if exact else [],
callback=callback)
errMsg = gdal.GetLastErrorMsg()
if errMsg == 'User terminated':
raise TimeoutError('raster2polygon timed out!')
raise TimeoutError('raster2polygon timed out!') if PY3 else TimeoutError_comp('raster2polygon timed out!')
if result is None:
raise Exception(errMsg)
# extract polygon
mem_layer.SetAttributeFilter('DN = 1')
element = mem_layer.GetNextFeature() # asserts that
shplyPoly = loads(element.GetGeometryRef().ExportToWkb())
mem_layer.SetAttributeFilter('DN = %s' %DN2extract)
from geopandas import GeoDataFrame
featCount = mem_layer.GetFeatureCount()
if maxfeatCount and featCount > maxfeatCount:
raise RuntimeError('Found %s features with DN=%s but maximum feature count was set to %s.'
%(featCount, DN2extract, maxfeatCount))
#if featCount>1:
#tmp = np.full((featCount,2), DN, geoArr.dtype)
#tmp[:,0] = range(featCount)
#GDF = GeoDataFrame(tmp, columns=['idx','DN'])
#def get_shplyPoly(GDF_row):
# if not is_timed_out(3):
# element = mem_layer.GetNextFeature()
# shplyPoly = loads(element.GetGeometryRef().ExportToWkb()).buffer(0)
# element = None
# return shplyPoly
# else:
# raise TimeoutError
#GDF['geometry'] = GDF.apply(get_shplyPoly, axis=1)
GDF = GeoDataFrame(columns=['geometry', 'DN'])
for i in range(featCount):
if not is_timed_out(timeout):
element = mem_layer.GetNextFeature()
GDF.loc[i] = [loads(element.GetGeometryRef().ExportToWkb()).buffer(0), DN2extract]
element = None
else:
raise TimeoutError('raster2polygon timed out!') if PY3 else TimeoutError_comp('raster2polygon timed out!')
GDF = GDF.dissolve(by='DN')
mem_ds = mem_layer = None
shplyPoly = GDF.loc[1,'geometry']
return shplyPoly
......@@ -4,8 +4,10 @@ __author__ = "Daniel Scheffler"
import math
import warnings
import numpy as np
from shapely.geometry import shape, Polygon, box
from geopandas import GeoDataFrame
from shapely.geometry import shape, Polygon, box, Point
from ..coord_trafo import mapYX2imYX
from ..coord_grid import find_nearest_grid_coord
......@@ -85,4 +87,23 @@ def find_line_intersection_point(line1, line2):
d = (det(*line1), det(*line2))
x = det(d, xdiff) / div
y = det(d, ydiff) / div
return x, y
\ No newline at end of file
return x, y
def polyVertices_outside_poly(inner_poly, outer_poly):
"""Checks if a shapely polygon (inner_poly) contains vertices that do not intersect 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
"""
if inner_poly.within(outer_poly):
# all vertices are inside outer_poly
return False
elif inner_poly.intersects(outer_poly):
# check if all vertices intersect with outer poly
GDF = GeoDataFrame(np.swapaxes(np.array(inner_poly.exterior.coords.xy), 0, 1), columns=['X', 'Y'])
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
\ No newline at end of file
......@@ -6,9 +6,11 @@ import numpy as np
import os
import warnings
from matplotlib import pyplot as plt
from six import PY3
# custom
from shapely.geometry import Polygon, box
from shapely.wkt import loads as shply_loads
from osgeo import gdal_array
# mpl_toolkits.basemap -> imported when GeoArray.show_map() is used
# dill -> imported when dumping GeoArray
......@@ -26,12 +28,12 @@ from ...geo.coord_grid import snap_bounds_to_pixGrid
from ...geo.coord_trafo import mapXY2imXY, imXY2mapXY, transform_any_prj, reproject_shapelyGeometry
from ...geo.projection import prj_equal, WKT2EPSG, EPSG2WKT
from ...geo.raster.conversion import raster2polygon
from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon
from ...geo.vector.topology import get_overlap_polygon, get_footprint_polygon, polyVertices_outside_poly
from ...geo.vector.geometry import boxObj
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...numeric.array import find_noDataVal, get_outFillZeroSaturated
from ...compatibility.python.exceptions import TimeoutError
from ...compatibility.python.exceptions import TimeoutError as TimeoutError_comp, \
FileNotFoundError as FileNotFoundError_comp
def _alias_property(key):
......@@ -86,10 +88,7 @@ class GeoArray(object):
assert ' ' not in self.arg, "The given path contains whitespaces. This is not supported by GDAL."
if not os.path.exists(self.filePath):
try:
raise FileNotFoundError(self.filePath)
except NameError: # FileNotFoundError is not available in Python 2.7
raise IOError(self.filePath)
raise FileNotFoundError(self.filePath) if PY3 else FileNotFoundError_comp(self.filePath)
ds = gdal.Open(self.filePath)
if not ds:
......@@ -289,17 +288,37 @@ class GeoArray(object):
self._footprint_poly = self.box.mapPoly
else:
try:
self._footprint_poly = raster2polygon(self, exact=False, progress=self.progress, q=self.q, timeout=10)
except TimeoutError:
self._footprint_poly = raster2polygon(self, exact=False, progress=self.progress, q=self.q,
maxfeatCount=10, timeout=3)
except (RuntimeError, TimeoutError, TimeoutError_comp):
if not self.q:
warnings.warn("\nCalculation of footprint polygon failed for %s '%s'. Using outer bounds. One "
"reason could be that the nodata value appears in the middle of the actual image. "
"reason could be that the nodata value appears within of the actual image. "
"To avoid this use another nodata value. Current nodata value is %s."
%(self.__class__.__name__, self.basename, self.nodata))
self._footprint_poly = self.box.mapPoly
# validation
assert not polyVertices_outside_poly(self._footprint_poly, self.box.mapPoly)
#assert self._footprint_poly
# for XY in self.corner_coord:
# assert self.GeoArray.box.mapPoly.contains(Point(XY)) or self.GeoArray.box.mapPoly.touches(Point(XY)), \
# "The corner position '%s' is outside of the %s." % (XY, self.imName)
return self._footprint_poly
@footprint_poly.setter
def footprint_poly(self, poly):
if isinstance(poly, Polygon):
self._footprint_poly = poly
elif isinstance(poly, str):
self._footprint_poly = shply_loads(poly)
else:
raise ValueError("'footprint_poly' can only be set from a shapely polygon or a WKT string.")
def __getitem__(self, given):
# TODO check if array cache contains the needed slice and return data from there
......
......@@ -9,13 +9,16 @@ _time_start = None
def is_timed_out(timeout):
global _time_start
if _time_start is None and timeout:
_time_start = time()
if timeout and time() - _time_start >= timeout:
_time_start = None
return True
if timeout is not None:
global _time_start
if _time_start is None:
_time_start = time()
if time() - _time_start >= timeout:
_time_start = None
return True
else:
return False
else:
return False
......
Supports Markdown
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