simplifier.py

import argparse
import sys
import os
from osgeo import ogr
from osgeo import osr
import anyjson
import shapely.geometry
import shapely.ops
import codecs
import time


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

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


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

in_field_count = in_defn.GetFieldCount()
for fld_index in range(in_field_count):
    src_fd = in_defn.GetFieldDefn( fld_index )
    fd = ogr.FieldDefn( src_fd.GetName(), src_fd.GetType() )
    fd.SetWidth( src_fd.GetWidth() )
    fd.SetPrecision( src_fd.GetPrecision() )
    shp_layer.CreateField( fd )


# Load geometries
geometries = []
for feature in in_layer:
  geometry = feature.GetGeometryRef()
  geometryType = geometry.GetGeometryType()

  if geometryType == ogr.wkbPolygon or geometryType == ogr.wkbMultiPolygon:
    shapelyGeometry = shapely.wkb.loads( geometry.ExportToWkb() )
    #if not shapelyGeometry.is_valid:
      #buffer to fix selfcrosses
      #shapelyGeometry = shapelyGeometry.buffer(0)
    if shapelyGeometry:
      geometries.append(shapelyGeometry)
in_layer.ResetReading()

start = int(round(time.time() * 1000))
# Simplification
points = []
connections = {}
counter = 0
for geom in geometries:
  counter += 1
  polygons = []

  if isinstance(geom, shapely.geometry.Polygon):
    polygons.append(geom)
  else:
    for polygon in geom:
      polygons.append(polygon)

  for polygon in polygons:
    if polygon.area > 0:
      lines = []
      lines.append(polygon.exterior)
      for line in polygon.interiors:
        lines.append(line)

      for line in lines:
        for i in range(len(line.coords)-1):
          indexFrom = i
          indexTo = i+1
          pointFrom = format % line.coords[indexFrom]
          pointTo = format % line.coords[indexTo]
          if pointFrom == pointTo:
            continue
          if not (pointFrom in connections):
            connections[pointFrom] = {}
          connections[pointFrom][pointTo] = 1
          if not (pointTo in connections):
            connections[pointTo] = {}
          connections[pointTo][pointFrom] = 1
print int(round(time.time() * 1000)) - start

simplifiedLines = {}
pivotPoints = {}
def simplifyRing(ring):
  coords = list(ring.coords)[0:-1]
  simpleCoords = []

  isPivot = False
  pointIndex = 0
  while not isPivot and pointIndex < len(coords):
    pointStr = format % coords[pointIndex]
    pointIndex += 1
    isPivot = ((len(connections[pointStr]) > 2) or (pointStr in pivotPoints))
  pointIndex = pointIndex - 1

  if not isPivot:
    simpleRing = shapely.geometry.LineString(coords).simplify(tolerance)
    if len(simpleRing.coords) <= 2:
      return None
    else:
      pivotPoints[format % coords[0]] = True
      pivotPoints[format % coords[-1]] = True
      simpleLineKey = format % coords[0]+':'+format % coords[1]+':'+format % coords[-1]
      simplifiedLines[simpleLineKey] = simpleRing.coords
      return simpleRing
  else:
    points = coords[pointIndex:len(coords)]
    points.extend(coords[0:pointIndex+1])
    iFrom = 0
    for i in range(1, len(points)):
      pointStr = format % points[i]
      if ((len(connections[pointStr]) > 2) or (pointStr in pivotPoints)):
        line = points[iFrom:i+1]
        lineKey = format % line[-1]+':'+format % line[-2]+':'+format % line[0]
        if lineKey in simplifiedLines:
          simpleLine = simplifiedLines[lineKey]
          simpleLine = list(reversed(simpleLine))
        else:
          simpleLine = shapely.geometry.LineString(line).simplify(tolerance).coords
          lineKey = format % line[0]+':'+format % line[1]+':'+format % line[-1]
          simplifiedLines[lineKey] = simpleLine
        simpleCoords.extend( simpleLine[0:-1] )
        iFrom = i
    if len(simpleCoords) <= 2:
      return None
    else:
      return shapely.geometry.LineString(simpleCoords)


def simplifyPolygon(polygon):
  simpleExtRing = simplifyRing(polygon.exterior)
  if simpleExtRing is None:
    return None
  simpleIntRings = []
  for ring in polygon.interiors:
    simpleIntRing = simplifyRing(ring)
    if simpleIntRing is not None:
      simpleIntRings.append(simpleIntRing)
  return shapely.geometry.Polygon(simpleExtRing, simpleIntRings)


results = []
for geom in geometries:
  polygons = []
  simplePolygons = []

  if isinstance(geom, shapely.geometry.Polygon):
    polygons.append(geom)
  else:
    for polygon in geom:
      polygons.append(polygon)

  for polygon in polygons:
    simplePolygon = simplifyPolygon(polygon)
    if not (simplePolygon is None or simplePolygon._geom is None):
      simplePolygons.append(simplePolygon)

  if len(simplePolygons) > 0:
    results.append(shapely.geometry.MultiPolygon(simplePolygons))
  else:
    results.append(None)


# Process all features in input layer.
in_feat = in_layer.GetNextFeature()
counter = 0
while in_feat is not None:
  if results[counter] is not None:
    out_feat = ogr.Feature( feature_def = shp_layer.GetLayerDefn() )
    out_feat.SetFrom( in_feat )
    out_feat.SetGeometryDirectly(
      ogr.CreateGeometryFromWkb(
        shapely.wkb.dumps(
          results[counter]
        )
      )
    )
    shp_layer.CreateFeature( out_feat )
    out_feat.Destroy()
  else:
    print 'geometry is too small: '+in_feat.GetField(16)

  in_feat.Destroy()
  in_feat = in_layer.GetNextFeature()
  counter += 1


# Cleanup
shp_ds.Destroy()
in_ds.Destroy()

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