15
infile = '/Users/kirilllebedev/Maps/50m-admin-0-countries/ne_50m_admin_0_countries.shp'
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()
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() )
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 )
44
for feature in in_layer:
45
geometry = feature.GetGeometryRef()
46
geometryType = geometry.GetGeometryType()
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)
54
geometries.append(shapelyGeometry)
55
in_layer.ResetReading()
57
start = int(round(time.time() * 1000))
62
for geom in geometries:
66
if isinstance(geom, shapely.geometry.Polygon):
70
polygons.append(polygon)
72
for polygon in polygons:
75
lines.append(polygon.exterior)
76
for line in polygon.interiors:
80
for i in range(len(line.coords)-1):
83
pointFrom = format % line.coords[indexFrom]
84
pointTo = format % line.coords[indexTo]
85
if pointFrom == pointTo:
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
97
def simplifyRing(ring):
98
coords = list(ring.coords)[0:-1]
103
while not isPivot and pointIndex < len(coords):
104
pointStr = format % coords[pointIndex]
106
isPivot = ((len(connections[pointStr]) > 2) or (pointStr in pivotPoints))
107
pointIndex = pointIndex - 1
110
simpleRing = shapely.geometry.LineString(coords).simplify(tolerance)
111
if len(simpleRing.coords) <= 2:
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
120
points = coords[pointIndex:len(coords)]
121
points.extend(coords[0:pointIndex+1])
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))
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] )
137
if len(simpleCoords) <= 2:
140
return shapely.geometry.LineString(simpleCoords)
143
def simplifyPolygon(polygon):
144
simpleExtRing = simplifyRing(polygon.exterior)
145
if simpleExtRing is None:
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)
156
for geom in geometries:
160
if isinstance(geom, shapely.geometry.Polygon):
161
polygons.append(geom)
164
polygons.append(polygon)
166
for polygon in polygons:
167
simplePolygon = simplifyPolygon(polygon)
168
if not (simplePolygon is None or simplePolygon._geom is None):
169
simplePolygons.append(simplePolygon)
171
if len(simplePolygons) > 0:
172
results.append(shapely.geometry.MultiPolygon(simplePolygons))
177
# Process all features in input layer.
178
in_feat = in_layer.GetNextFeature()
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(
191
shp_layer.CreateFeature( out_feat )
194
print 'geometry is too small: '+in_feat.GetField(16)
197
in_feat = in_layer.GetNextFeature()
205
print int(round(time.time() * 1000)) - start