conversion.py 6.45 KB
Newer Older
1
# -*- coding: utf-8 -*-
2

3
4
# 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.
5
#
6
7
8
9
# 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/)
10
11
12
13
14
#
# 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).
#
15
16
17
# 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
18
#
19
#   http://www.apache.org/licenses/LICENSE-2.0
20
#
21
22
23
24
25
# 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.
26

27
28
29
__author__ = "Daniel Scheffler"

from shapely.wkb import loads
30
from osgeo import gdal, osr, ogr
31

32
33
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...processing.progress_mon import ProgressBar, Timer
34
from ..raster.reproject import warp_ndarray
35
from ..vector.topology import get_bounds_polygon, polyVertices_outside_poly, get_overlap_polygon
36
37


38
def raster2polygon(array, gt, prj, DN2extract=1, exact=True, maxfeatCount=None,
Daniel Scheffler's avatar
Daniel Scheffler committed
39
                   timeout=None, progress=True, q=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
40
    """Calculate a footprint polygon for the given array.
41

42
    :param array:             2D numpy array
43
44
    :param gt:                GDAL GeoTransform
    :param prj:               projection as WKT string, 'EPSG:1234' or <EPSG_int>
45
    :param DN2extract:        <int, float> pixel value to create polygons for
46
    :param exact:             whether to compute the exact footprint polygon or a simplified one for speed
47
48
    :param maxfeatCount:      <int> the maximum expected number of polygons. If more polygons are found, every further
                              processing is cancelled and a RunTimeError is raised.
49
50
51
    :param timeout:           breaks the process after a given time in seconds
    :param progress:          show progress bars (default: True)
    :param q:                 quiet mode (default: False)
52
53
    :return:
    """
54
    assert array.ndim == 2, "Only 2D arrays are supported. Got a %sD array." % array.ndim
55
56
    gt_orig = gt
    shape_orig = array.shape
57

58
    # downsample input array in case it has more than 1e8 pixels to prevent crash
59
    if not exact and array.size > 1e8:  # 10000 x 10000 px
60
        # downsample with nearest neighbour
61
62
        zoom_factor = (8000 * 8000 / array.size) ** 0.5
        array, gt, prj = warp_ndarray(array, gt, prj,
Daniel Scheffler's avatar
Daniel Scheffler committed
63
64
                                      out_gsd=(gt[1] / zoom_factor,
                                               gt[5] / zoom_factor),
65
66
                                      rspAlg='near',
                                      q=True)
67

68
    src_ds = get_GDAL_ds_inmem(array, gt, prj)
69
70
71
72
    src_band = src_ds.GetRasterBand(1)

    # Create a memory OGR datasource to put results in.
    mem_drv = ogr.GetDriverByName('Memory')
73
    mem_ds = mem_drv.CreateDataSource('out')
74
75

    srs = osr.SpatialReference()
76
    srs.ImportFromWkt(prj)
77
78
79
80
81
82

    mem_layer = mem_ds.CreateLayer('poly', srs, ogr.wkbPolygon)

    fd = ogr.FieldDefn('DN', ogr.OFTInteger)
    mem_layer.CreateField(fd)

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    # set callback
    callback = \
        ProgressBar(prefix='Polygonize progress    ',
                    suffix='Complete',
                    barLength=50,
                    timeout=timeout,
                    use_as_callback=True) \
        if progress and not q else Timer(timeout, use_as_callback=True) if timeout else None

    # run the algorithm
    status = gdal.Polygonize(src_band,
                             src_band.GetMaskBand(),
                             mem_layer,
                             0,
                             ["8CONNECTED=8"] if exact else [],
                             callback=callback)

    # handle exit status other than 0 (fail)
    if status != 0:
        errMsg = gdal.GetLastErrorMsg()

        # Catch the KeyboardInterrupt raised in case of a timeout within the callback. It seems like there is no other
        # way to catch exceptions within callbacks.
        if errMsg == 'User terminated':
            raise TimeoutError('raster2polygon timed out!')

        raise Exception(errMsg)
Daniel Scheffler's avatar
Daniel Scheffler committed
110

111
    # extract polygon
112
    mem_layer.SetAttributeFilter('DN = %s' % DN2extract)
Daniel Scheffler's avatar
Daniel Scheffler committed
113
114
115
116

    from geopandas import GeoDataFrame
    featCount = mem_layer.GetFeatureCount()

117
    if not featCount:
118
        raise RuntimeError('No features with DN=%s found in the input image.' % DN2extract)
Daniel Scheffler's avatar
Daniel Scheffler committed
119
120
    if maxfeatCount and featCount > maxfeatCount:
        raise RuntimeError('Found %s features with DN=%s but maximum feature count was set to %s.'
121
                           % (featCount, DN2extract, maxfeatCount))
Daniel Scheffler's avatar
Daniel Scheffler committed
122

123
124
125
    # tmp = np.full((featCount,2), DN, geoArr.dtype)
    # tmp[:,0] = range(featCount)
    # GDF = GeoDataFrame(tmp, columns=['idx','DN'])
Daniel Scheffler's avatar
Daniel Scheffler committed
126

127
    # def get_shplyPoly(GDF_row):
Daniel Scheffler's avatar
Daniel Scheffler committed
128
129
130
131
132
133
134
135
    #    if not is_timed_out(3):
    #        element   = mem_layer.GetNextFeature()
    #        shplyPoly = loads(element.GetGeometryRef().ExportToWkb()).buffer(0)
    #        element   = None
    #        return shplyPoly
    #    else:
    #        raise TimeoutError

136
    # GDF['geometry'] = GDF.apply(get_shplyPoly, axis=1)
Daniel Scheffler's avatar
Daniel Scheffler committed
137
138

    GDF = GeoDataFrame(columns=['geometry', 'DN'])
139
    GDF.DN = GDF.DN.astype(float)
140
    timer = Timer(timeout)
Daniel Scheffler's avatar
Daniel Scheffler committed
141
    for i in range(featCount):
142
        if not timer.timed_out:
143
            element = mem_layer.GetNextFeature()
144
145
            wkb = bytes(element.GetGeometryRef().ExportToWkb())
            GDF.loc[i] = [loads(wkb).buffer(0), DN2extract]
146
            del element
Daniel Scheffler's avatar
Daniel Scheffler committed
147
        else:
148
            raise TimeoutError('raster2polygon timed out!')
Daniel Scheffler's avatar
Daniel Scheffler committed
149
150
151

    GDF = GDF.dissolve(by='DN')

152
    del mem_ds, mem_layer
153

154
    shplyPoly = GDF.loc[1, 'geometry']
155
156
157
158
159
160
161
162
163

    # 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']

164
    return shplyPoly