conversion.py 5.33 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# py_tools_ds
#
# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@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).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# 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/>.

24
25
__author__ = "Daniel Scheffler"

Daniel Scheffler's avatar
Daniel Scheffler committed
26
from six import PY3
27
28

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

31
32
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...processing.progress_mon import ProgressBar, Timer
Daniel Scheffler's avatar
Daniel Scheffler committed
33
from ...compatibility.python.exceptions import TimeoutError as TimeoutError_comp
34
35


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

40
    :param array:             2D numpy array
41
42
    :param gt:
    :param prj:
43
    :param DN2extract:        <int, float> pixel value to create polygons for
44
    :param exact:
45
46
    :param maxfeatCount:      <int> the maximum expected number of polygons. If more polygons are found, every further
                              processing is cancelled and a RunTimeError is raised.
47
48
49
    :param timeout:           breaks the process after a given time in seconds
    :param progress:          show progress bars (default: True)
    :param q:                 quiet mode (default: False)
50
51
    :return:
    """
52
53
    assert array.ndim == 2, "Only 2D arrays are supported. Got a %sD array." % array.ndim

54
55
    # downsample input array in case is has more than 1e8 pixels to prevent crash
    if not exact and array.size > 1e8:  # 10000 x 10000 px
56
        zoom_factor = 1e8 / array.size
57

58
        # downsample with nearest neighbour
59
        from skimage.transform import rescale  # import here to avoid static TLS import error
60
        array = rescale(array, zoom_factor, order=0, preserve_range=True, mode='edge').astype(bool)
61
62
63
64
65
66

        # update pixel sizes within gt
        gt = list(gt)
        gt[1] /= zoom_factor
        gt[5] /= zoom_factor

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

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

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

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

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

82
    # set callback
Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
83
84
    callback = ProgressBar(prefix='Polygonize progress    ', suffix='Complete', barLength=50, timeout=timeout,
                           use_as_callback=True) \
85
        if progress and not q else Timer(timeout, use_as_callback=True) if timeout else None
86
87

    # run the algorithm
Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
88
    status = gdal.Polygonize(src_band, src_band.GetMaskBand(), mem_layer, 0, ["8CONNECTED=8"] if exact else [],
89
                             callback=callback)
Daniel Scheffler's avatar
Daniel Scheffler committed
90

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
91
92
93
94
95
    # handle exit status other than 0 (fail)
    if status != 0:
        errMsg = gdal.GetLastErrorMsg()
        if errMsg == 'User terminated':
            raise TimeoutError('raster2polygon timed out!') if PY3 else TimeoutError_comp('raster2polygon timed out!')
96
        raise Exception(errMsg)
97
98

    # extract polygon
99
    mem_layer.SetAttributeFilter('DN = %s' % DN2extract)
Daniel Scheffler's avatar
Daniel Scheffler committed
100
101
102
103

    from geopandas import GeoDataFrame
    featCount = mem_layer.GetFeatureCount()

104
    if not featCount:
105
        raise RuntimeError('No features with DN=%s found in the input image.' % DN2extract)
Daniel Scheffler's avatar
Daniel Scheffler committed
106
107
    if maxfeatCount and featCount > maxfeatCount:
        raise RuntimeError('Found %s features with DN=%s but maximum feature count was set to %s.'
108
                           % (featCount, DN2extract, maxfeatCount))
Daniel Scheffler's avatar
Daniel Scheffler committed
109

110
111
112
    # 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
113

114
    # def get_shplyPoly(GDF_row):
Daniel Scheffler's avatar
Daniel Scheffler committed
115
116
117
118
119
120
121
122
    #    if not is_timed_out(3):
    #        element   = mem_layer.GetNextFeature()
    #        shplyPoly = loads(element.GetGeometryRef().ExportToWkb()).buffer(0)
    #        element   = None
    #        return shplyPoly
    #    else:
    #        raise TimeoutError

123
    # GDF['geometry'] = GDF.apply(get_shplyPoly, axis=1)
Daniel Scheffler's avatar
Daniel Scheffler committed
124
125

    GDF = GeoDataFrame(columns=['geometry', 'DN'])
126
    timer = Timer(timeout)
Daniel Scheffler's avatar
Daniel Scheffler committed
127
    for i in range(featCount):
128
        if not timer.timed_out:
129
            element = mem_layer.GetNextFeature()
Daniel Scheffler's avatar
Daniel Scheffler committed
130
            GDF.loc[i] = [loads(element.GetGeometryRef().ExportToWkb()).buffer(0), DN2extract]
131
            del element
Daniel Scheffler's avatar
Daniel Scheffler committed
132
133
134
135
136
        else:
            raise TimeoutError('raster2polygon timed out!') if PY3 else TimeoutError_comp('raster2polygon timed out!')

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

137
    del mem_ds, mem_layer
138

139
    shplyPoly = GDF.loc[1, 'geometry']
140
    return shplyPoly