conversion.py 5.47 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
import numpy as np
30

31
32
33
34
35
36
37
38
39
try:
    import gdal
    import ogr
    import osr
except ImportError:
    from osgeo import gdal
    from osgeo import ogr
    from osgeo import osr

40
41
from ...io.raster.gdal import get_GDAL_ds_inmem
from ...processing.progress_mon import ProgressBar, Timer
Daniel Scheffler's avatar
Daniel Scheffler committed
42
from ...compatibility.python.exceptions import TimeoutError as TimeoutError_comp
43
44


45
def raster2polygon(array, gt, prj, DN2extract=1, exact=True, maxfeatCount=None,
Daniel Scheffler's avatar
Daniel Scheffler committed
46
                   timeout=None, progress=True, q=False):
47
    """Calculates a footprint polygon for the given array.
48

49
    :param array:             2D numpy array
50
51
    :param gt:
    :param prj:
52
    :param DN2extract:        <int, float> pixel value to create polygons for
53
    :param exact:
54
55
    :param maxfeatCount:      <int> the maximum expected number of polygons. If more polygons are found, every further
                              processing is cancelled and a RunTimeError is raised.
56
57
58
    :param timeout:           breaks the process after a given time in seconds
    :param progress:          show progress bars (default: True)
    :param q:                 quiet mode (default: False)
59
60
61
    :return:
    """

62
63
    assert array.ndim == 2, "Only 2D arrays are supported. Got a %sD array." % array.ndim

64
65
    # 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
66
67
68
        zoom_factor = 0.5

        # downsample to half size, nearest neighbour
69
        from skimage.transform import rescale  # import here to avoid static TLS import error
70
71
72
73
74
75
76
        array = rescale(array, zoom_factor, order=0, preserve_range=True, mode='edge').astype(np.bool)

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

77
    src_ds = get_GDAL_ds_inmem(array, gt, prj)
78
79
80
81
    src_band = src_ds.GetRasterBand(1)

    # Create a memory OGR datasource to put results in.
    mem_drv = ogr.GetDriverByName('Memory')
82
    mem_ds = mem_drv.CreateDataSource('out')
83
84

    srs = osr.SpatialReference()
85
    srs.ImportFromWkt(prj)
86
87
88
89
90
91

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

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

92
    # set callback
Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
93
94
    callback = ProgressBar(prefix='Polygonize progress    ', suffix='Complete', barLength=50, timeout=timeout,
                           use_as_callback=True) \
95
        if progress and not q else Timer(timeout, use_as_callback=True) if timeout else None
96
97

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

Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
101
102
103
104
105
    # 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!')
106
        raise Exception(errMsg)
107
108

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

    from geopandas import GeoDataFrame
    featCount = mem_layer.GetFeatureCount()

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

120
121
122
    # 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
123

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

133
    # GDF['geometry'] = GDF.apply(get_shplyPoly, axis=1)
Daniel Scheffler's avatar
Daniel Scheffler committed
134
135

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

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

147
    del mem_ds, mem_layer
148

149
    shplyPoly = GDF.loc[1, 'geometry']
150
    return shplyPoly