lavkach3

1
import argparse
2
import sys
3
import os
4
from osgeo import ogr
5
from osgeo import osr
6
import anyjson
7
import shapely.geometry
8
import shapely.ops
9
import codecs
10
import time
11

12

13
format = '%.8f %.8f'
14
tolerance = 0.01
15
infile = '/Users/kirilllebedev/Maps/50m-admin-0-countries/ne_50m_admin_0_countries.shp'
16
outfile = 'map.shp'
17

18
# Open the datasource to operate on.
19
in_ds = ogr.Open( infile, update = 0 )
20
in_layer = in_ds.GetLayer( 0 )
21
in_defn = in_layer.GetLayerDefn()
22

23

24
# Create output file with similar information.
25
shp_driver = ogr.GetDriverByName( 'ESRI Shapefile' )
26
if os.path.exists('map.shp'):
27
  shp_driver.DeleteDataSource( outfile )
28
shp_ds = shp_driver.CreateDataSource( outfile )
29
shp_layer = shp_ds.CreateLayer( in_defn.GetName(),
30
                                geom_type = in_defn.GetGeomType(),
31
                                srs = in_layer.GetSpatialRef() )
32

33
in_field_count = in_defn.GetFieldCount()
34
for fld_index in range(in_field_count):
35
    src_fd = in_defn.GetFieldDefn( fld_index )
36
    fd = ogr.FieldDefn( src_fd.GetName(), src_fd.GetType() )
37
    fd.SetWidth( src_fd.GetWidth() )
38
    fd.SetPrecision( src_fd.GetPrecision() )
39
    shp_layer.CreateField( fd )
40

41

42
# Load geometries
43
geometries = []
44
for feature in in_layer:
45
  geometry = feature.GetGeometryRef()
46
  geometryType = geometry.GetGeometryType()
47

48
  if geometryType == ogr.wkbPolygon or geometryType == ogr.wkbMultiPolygon:
49
    shapelyGeometry = shapely.wkb.loads( geometry.ExportToWkb() )
50
    #if not shapelyGeometry.is_valid:
51
      #buffer to fix selfcrosses
52
      #shapelyGeometry = shapelyGeometry.buffer(0)
53
    if shapelyGeometry:
54
      geometries.append(shapelyGeometry)
55
in_layer.ResetReading()
56

57
start = int(round(time.time() * 1000))
58
# Simplification
59
points = []
60
connections = {}
61
counter = 0
62
for geom in geometries:
63
  counter += 1
64
  polygons = []
65

66
  if isinstance(geom, shapely.geometry.Polygon):
67
    polygons.append(geom)
68
  else:
69
    for polygon in geom:
70
      polygons.append(polygon)
71

72
  for polygon in polygons:
73
    if polygon.area > 0:
74
      lines = []
75
      lines.append(polygon.exterior)
76
      for line in polygon.interiors:
77
        lines.append(line)
78

79
      for line in lines:
80
        for i in range(len(line.coords)-1):
81
          indexFrom = i
82
          indexTo = i+1
83
          pointFrom = format % line.coords[indexFrom]
84
          pointTo = format % line.coords[indexTo]
85
          if pointFrom == pointTo:
86
            continue
87
          if not (pointFrom in connections):
88
            connections[pointFrom] = {}
89
          connections[pointFrom][pointTo] = 1
90
          if not (pointTo in connections):
91
            connections[pointTo] = {}
92
          connections[pointTo][pointFrom] = 1
93
print int(round(time.time() * 1000)) - start
94

95
simplifiedLines = {}
96
pivotPoints = {}
97
def simplifyRing(ring):
98
  coords = list(ring.coords)[0:-1]
99
  simpleCoords = []
100

101
  isPivot = False
102
  pointIndex = 0
103
  while not isPivot and pointIndex < len(coords):
104
    pointStr = format % coords[pointIndex]
105
    pointIndex += 1
106
    isPivot = ((len(connections[pointStr]) > 2) or (pointStr in pivotPoints))
107
  pointIndex = pointIndex - 1
108

109
  if not isPivot:
110
    simpleRing = shapely.geometry.LineString(coords).simplify(tolerance)
111
    if len(simpleRing.coords) <= 2:
112
      return None
113
    else:
114
      pivotPoints[format % coords[0]] = True
115
      pivotPoints[format % coords[-1]] = True
116
      simpleLineKey = format % coords[0]+':'+format % coords[1]+':'+format % coords[-1]
117
      simplifiedLines[simpleLineKey] = simpleRing.coords
118
      return simpleRing
119
  else:
120
    points = coords[pointIndex:len(coords)]
121
    points.extend(coords[0:pointIndex+1])
122
    iFrom = 0
123
    for i in range(1, len(points)):
124
      pointStr = format % points[i]
125
      if ((len(connections[pointStr]) > 2) or (pointStr in pivotPoints)):
126
        line = points[iFrom:i+1]
127
        lineKey = format % line[-1]+':'+format % line[-2]+':'+format % line[0]
128
        if lineKey in simplifiedLines:
129
          simpleLine = simplifiedLines[lineKey]
130
          simpleLine = list(reversed(simpleLine))
131
        else:
132
          simpleLine = shapely.geometry.LineString(line).simplify(tolerance).coords
133
          lineKey = format % line[0]+':'+format % line[1]+':'+format % line[-1]
134
          simplifiedLines[lineKey] = simpleLine
135
        simpleCoords.extend( simpleLine[0:-1] )
136
        iFrom = i
137
    if len(simpleCoords) <= 2:
138
      return None
139
    else:
140
      return shapely.geometry.LineString(simpleCoords)
141

142

143
def simplifyPolygon(polygon):
144
  simpleExtRing = simplifyRing(polygon.exterior)
145
  if simpleExtRing is None:
146
    return None
147
  simpleIntRings = []
148
  for ring in polygon.interiors:
149
    simpleIntRing = simplifyRing(ring)
150
    if simpleIntRing is not None:
151
      simpleIntRings.append(simpleIntRing)
152
  return shapely.geometry.Polygon(simpleExtRing, simpleIntRings)
153

154

155
results = []
156
for geom in geometries:
157
  polygons = []
158
  simplePolygons = []
159

160
  if isinstance(geom, shapely.geometry.Polygon):
161
    polygons.append(geom)
162
  else:
163
    for polygon in geom:
164
      polygons.append(polygon)
165

166
  for polygon in polygons:
167
    simplePolygon = simplifyPolygon(polygon)
168
    if not (simplePolygon is None or simplePolygon._geom is None):
169
      simplePolygons.append(simplePolygon)
170

171
  if len(simplePolygons) > 0:
172
    results.append(shapely.geometry.MultiPolygon(simplePolygons))
173
  else:
174
    results.append(None)
175

176

177
# Process all features in input layer.
178
in_feat = in_layer.GetNextFeature()
179
counter = 0
180
while in_feat is not None:
181
  if results[counter] is not None:
182
    out_feat = ogr.Feature( feature_def = shp_layer.GetLayerDefn() )
183
    out_feat.SetFrom( in_feat )
184
    out_feat.SetGeometryDirectly(
185
      ogr.CreateGeometryFromWkb(
186
        shapely.wkb.dumps(
187
          results[counter]
188
        )
189
      )
190
    )
191
    shp_layer.CreateFeature( out_feat )
192
    out_feat.Destroy()
193
  else:
194
    print 'geometry is too small: '+in_feat.GetField(16)
195

196
  in_feat.Destroy()
197
  in_feat = in_layer.GetNextFeature()
198
  counter += 1
199

200

201
# Cleanup
202
shp_ds.Destroy()
203
in_ds.Destroy()
204

205
print int(round(time.time() * 1000)) - start