first mockup in go

static/bower_components/jvectormap/converter/simplifier.py

@@ -0,0 +1,205 @@
+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