File: geom.py

package info (click to toggle)
imposm 2.6.0+ds-6
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 996 kB
  • sloc: python: 4,164; cpp: 878; makefile: 13; sql: 1
file content (434 lines) | stat: -rw-r--r-- 15,222 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
# Copyright 2011, 2012 Omniscale (http://omniscale.com)
#
# 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
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# 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.

from __future__ import division

import math
import codecs
import os
import shapely.geometry
import shapely.geos
import shapely.prepared
from shapely.geometry.base import BaseGeometry
from shapely import geometry
from shapely import wkt
from shapely.ops import cascaded_union, linemerge
from shapely.topology import TopologicalError

try:
    import rtree
except ImportError:
    rtree = None

from imposm import config
from imposm.util.geom import load_polygons, load_datasource, build_multipolygon

import logging
log = logging.getLogger(__name__)

class InvalidGeometryError(Exception):
    pass

class IncompletePolygonError(Exception):
    pass


TOLERANCE_DEEGREES = 1e-8
TOLERANCE_METERS = 1e-3

# older versions had unhandled floating point execptions in .buffer(0)
#SHAPELY_SUPPORTS_BUFFER = shapely.geos.geos_capi_version >= (1, 6, 0)
# The geos_capi_version check fails with Shapely 1.5.12 & GEOS 3.5.0,
# GEOS >= 3.1.1 has CAPI >= 1.6.0 and that's been in Debian for a while.
SHAPELY_SUPPORTS_BUFFER = True

def validate_and_simplify(geom, meter_units=False):
    if SHAPELY_SUPPORTS_BUFFER:
        try:
            # buffer(0) is nearly fast as is_valid
            return geom.buffer(0)
        except ValueError:
            # shapely raises ValueError if buffer(0) result is empty
            raise InvalidGeometryError('geometry is empty')

    orig_geom = geom
    if not geom.is_valid:
        tolerance = TOLERANCE_METERS if meter_units else TOLERANCE_DEEGREES
        try:
            geom = geom.simplify(tolerance, False)
        except ValueError:
            # shapely raises ValueError if buffer(0) result is empty
            raise InvalidGeometryError('geometry is empty')

        if not geom.is_valid:
            raise InvalidGeometryError('geometry is invalid, could not simplify: %s' %
                                       orig_geom)
    return geom

class GeomBuilder(object):
    def build(self, osm_elem):
        # TODO is this method still in use?
        try:
            if isinstance(osm_elem.coords, BaseGeometry):
                return osm_elem.coords
            geom_wkt = self.to_wkt(osm_elem.coords)
            if geom_wkt is not None:
                geom = wkt.loads(geom_wkt)
        except Exception, ex:
            raise InvalidGeometryError('unable to build geometry %s: %s %s' %
                                       (osm_elem.osm_id, ex, osm_elem.coords))
        if geom_wkt is None or geom is None:
            # unable to build valid wkt (non closed polygon, etc)
            raise InvalidGeometryError()
        return geom

    def check_geom_type(self, geom):
        return

    def build_geom(self, osm_elem):
        try:
            if isinstance(osm_elem.coords, BaseGeometry):
                if osm_elem.coords.is_empty:
                    raise InvalidGeometryError('empty geometry')
                self.check_geom_type(osm_elem.coords)
                return osm_elem.coords
            geom = self.to_geom(osm_elem.coords)
        except InvalidGeometryError, ex:
            raise InvalidGeometryError('unable to build geometry %s: %s %s' %
                                       (osm_elem.osm_id, ex, osm_elem.coords))
        if geom is None:
            # unable to build valid wkt (non closed polygon, etc)
            raise InvalidGeometryError()
        return geom

class PointBuilder(GeomBuilder):
    def to_wkt(self, data):
        if len(data) != 2:
            return None
        return 'POINT(%f %f)' % data

    def to_geom(self, data):
        if len(data) != 2:
            return None
        return geometry.Point(*data)

    def check_geom_type(self, geom):
        if geom.type != 'Point':
            raise InvalidGeometryError('expected Point, got %s' % geom.type)

    def build_checked_geom(self, osm_elem, validate=False):
        geom = self.build_geom(osm_elem)
        if not validate or geom.is_valid:
            return geom
        else:
            raise InvalidGeometryError('invalid geometry for %s: %s, %s' %
                                       (osm_elem.osm_id, geom, osm_elem.coords))

class PolygonBuilder(GeomBuilder):
    def to_wkt(self, data):
        if len(data) >= 4 and data[0] == data[-1]:
            return 'POLYGON((' + ', '.join('%f %f' % p for p in data) + '))'
        return None

    def to_geom(self, data):
        if len(data) >= 4 and data[0] == data[-1]:
            return geometry.Polygon(data)
        return None

    def check_geom_type(self, geom):
        if geom.type not in ('Polygon', 'MultiPolygon'):
            raise InvalidGeometryError('expected Polygon or MultiPolygon, got %s' % geom.type)

    def build_checked_geom(self, osm_elem, validate=False):
        geom = self.build_geom(osm_elem)
        if not validate:
            return geom
        try:
            return validate_and_simplify(geom)
        except InvalidGeometryError:
            raise InvalidGeometryError('invalid geometry for %s: %s, %s' %
                                       (osm_elem.osm_id, geom, osm_elem.coords))

class LineStringBuilder(GeomBuilder):
    def to_wkt(self, data):
        if len(data) <= 1:
            return None
        if len(data) == 2 and data[0] == data[1]:
            return None
        return 'LINESTRING(' + ', '.join('%f %f' % p for p in data) + ')'

    def check_geom_type(self, geom):
        if geom.type != 'LineString':
            raise InvalidGeometryError('expected LineString, got %s' % geom.type)

    def to_geom(self, data, max_length=None):
        if len(data) <= 1:
            return None
        if len(data) == 2 and data[0] == data[1]:
            return None
        if max_length is None:
            max_length = config.imposm_linestring_max_length
        if max_length and len(data) > max_length:
            chunks = math.ceil(len(data) / max_length)
            length = int(len(data) // chunks)
            lines = []
            for i in xrange(1, len(data), length):
                lines.append(geometry.LineString(data[i-1:i+length]))
            return lines
        return geometry.LineString(data)

    def build_checked_geom(self, osm_elem, validate=False):
        geom = self.build_geom(osm_elem)
        if not validate or geom.is_valid:
            return geom
        else:
            raise InvalidGeometryError('invalid geometry for %s: %s, %s' %
                                       (osm_elem.osm_id, geom, osm_elem.coords))

def tile_bbox(bbox, grid_width):
    """
    Tile bbox into multiple sub-boxes, each of `grid_width` size.

    >>> list(tile_bbox((-1, 1, 0.49, 1.51), 0.5)) #doctest: +NORMALIZE_WHITESPACE
    [(-1.0, 1.0, -0.5, 1.5),
     (-1.0, 1.5, -0.5, 2.0),
     (-0.5, 1.0, 0.0, 1.5),
     (-0.5, 1.5, 0.0, 2.0),
     (0.0, 1.0, 0.5, 1.5),
     (0.0, 1.5, 0.5, 2.0)]
    """
    min_x = math.floor(bbox[0]/grid_width) * grid_width
    min_y = math.floor(bbox[1]/grid_width) * grid_width
    max_x = math.ceil(bbox[2]/grid_width) * grid_width
    max_y = math.ceil(bbox[3]/grid_width) * grid_width

    x_steps = math.ceil((max_x - min_x) / grid_width)
    y_steps = math.ceil((max_y - min_y) / grid_width)

    for x in xrange(int(x_steps)):
        for y in xrange(int(y_steps)):
            yield (
                min_x + x * grid_width,
                min_y + y * grid_width,
                min_x + (x + 1) * grid_width,
                min_y + (y + 1)* grid_width,
            )

def split_polygon_at_grid(geom, grid_width=0.1, current_grid_width=10.0):
    """
    >>> p = list(split_polygon_at_grid(geometry.box(-0.5, 1, 0.2, 2), 1))
    >>> p[0].contains(geometry.box(-0.5, 1, 0, 2))
    True
    >>> p[0].area == geometry.box(-0.5, 1, 0, 2).area
    True
    >>> p[1].contains(geometry.box(0, 1, 0.2, 2))
    True
    >>> p[1].area == geometry.box(0, 1, 0.2, 2).area
    True
    """
    if not geom.is_valid:
        geom = geom.buffer(0)

    for i, split_box in enumerate(tile_bbox(geom.bounds, current_grid_width)):
        try:
            polygon_part = geom.intersection(shapely.geometry.box(*split_box))
        except TopologicalError:
            continue
        if not polygon_part.is_empty and polygon_part.type.endswith('Polygon'):
            if grid_width >= current_grid_width:
                yield polygon_part
            else:
                for part in split_polygon_at_grid(polygon_part, grid_width, current_grid_width/10.0):
                    yield part

def load_geom(source):
    geom = load_datasource(source)
    if geom:
        # get the first and maybe only geometry
        if not check_wgs84_srs(geom[0]):
            log.error('Geometry is not in EPSG:4326')
            return None
        if rtree:
            return LimitRTreeGeometry(geom)
        else:
            log.info('You should install RTree for large --limit-to polygons')
            return LimitPolygonGeometry(build_multipolygon(geom)[1])
    return None

def check_wgs84_srs(geom):
    bbox = geom.bounds
    if bbox[0] >= -180 and bbox[1] >= -90 and bbox[2] <= 180 and bbox[3] <= 90:
        return True
    return False

class EmtpyGeometryError(Exception):
    pass

class LimitPolygonGeometry(object):
    def __init__(self, shapely_geom):
        self._geom = shapely_geom
        self._prepared_geom = None
        self._prepared_counter = 0
        self._prepared_max = 100000

    @property
    def geom(self):
        # GEOS internal data structure for prepared geometries grows over time,
        # recreate to limit memory consumption
        if not self._prepared_geom or self._prepared_counter > self._prepared_max:
            self._prepared_geom = shapely.prepared.prep(self._geom)
            self._prepared_counter = 0
        self._prepared_counter += 1
        return self._prepared_geom

    def intersection(self, geom):
        if self.geom.contains_properly(geom):
            # no need to limit contained geometries
            return geom

        new_geom = None
        if self.geom.intersects(geom):
            try:
                # can not use intersection with prepared geom
                new_geom = self._geom.intersection(geom)
            except TopologicalError:
                pass

        if not new_geom or new_geom.is_empty:
            raise EmtpyGeometryError('No intersection or empty geometry')

        new_geom = filter_geometry_by_type(new_geom, geom.type)
        if new_geom:
            return new_geom

        raise EmtpyGeometryError('No intersection or empty geometry')

def filter_geometry_by_type(geometry, geom_type):
    """
    Filter (multi)geometry for compatible `geom_type`,
    because we can't insert points into linestring tables for example
    """
    if geometry.type == geom_type:
        # same type is fine
        return geometry

    if geometry.type == 'Polygon' and geom_type == 'MultiPolygon':
        # multipolygon mappings also support polygons
        return geometry

    if geometry.type == 'MultiPolygon' and geom_type == 'Polygon':
        # polygon mappings should also support multipolygons
        return geometry

    if hasattr(geometry, 'geoms'):
        # GeometryCollection or MultiLineString? return list of geometries
        geoms = []
        for part in geometry.geoms:
            # only parts with same type
            if part.type == geom_type:
                geoms.append(part)

        if geoms:
            return geoms

    return None

def flatten_polygons(polygons):
    for polygon in polygons:
        if polygon.type == 'MultiPolygon':
            for p in polygon.geoms:
                yield p
        else:
            yield polygon

def flatten_linestrings(linestrings):
    for linestring in linestrings:
        if linestring.type == 'MultiLineString':
            for ls in linestring.geoms:
                yield ls
        else:
            yield linestring

def filter_invalid_linestrings(linestrings):
    for linestring in linestrings:
        # filter out tiny linestrings, can become invalid geometries in postgis
        if linestring.length > 1e-9:
            yield linestring

class LimitRTreeGeometry(object):
    def __init__(self, polygons):
        index = rtree.index.Index()
        sub_polygons = []
        part_idx = 0
        for polygon in polygons:
            for part in split_polygon_at_grid(polygon):
                sub_polygons.append(part)
                index.insert(part_idx, part.bounds)
                part_idx += 1

        self.polygons = sub_polygons
        self.index = index

    def intersection(self, geom):
        intersection_ids = list(self.index.intersection(geom.bounds))

        if not intersection_ids:
            raise EmtpyGeometryError('No intersection or empty geometry')

        intersections = []
        for i in intersection_ids:

            polygon = self.polygons[i]

            if polygon.contains(geom):
                return geom

            if polygon.intersects(geom):
                try:
                    new_geom_part = polygon.intersection(geom)
                    new_geom_part = filter_geometry_by_type(new_geom_part, geom.type)
                    if new_geom_part:
                        if isinstance(new_geom_part, list):
                            intersections.extend(new_geom_part)
                        else:
                            intersections.append(new_geom_part)
                except TopologicalError:
                        pass

        if not intersections:
            raise EmtpyGeometryError('No intersection or empty geometry')

        # intersections from multiple sub-polygons
        # try to merge them back to a single geometry
        try:
            if geom.type.endswith('Polygon'):
                union = cascaded_union(list(flatten_polygons(intersections)))
            elif geom.type.endswith('LineString'):
                linestrings = flatten_linestrings(intersections)
                linestrings = list(filter_invalid_linestrings(linestrings))
                if not linestrings:
                    raise EmtpyGeometryError()
                union = linemerge(linestrings)
                if union.type == 'MultiLineString':
                    union = list(union.geoms)
            elif geom.type == 'Point':
                union = intersections[0]
            else:
                log.warn('unexpexted geometry type %s', geom.type)
                raise EmtpyGeometryError()
        except ValueError, ex:
            # likely an 'No Shapely geometry can be created from null value' error
            log.warn('could not create union: %s', ex)
            raise EmtpyGeometryError()
        return union