Merge pull request #21 from calvinmetcalf/topo

TopoJSON support
project-open-data · 2400d16efc6cfb36554cd64e07c90ba544b2099c · feomike committed Mar 24, 2015 · Mar 24, 2015 · 2400d16
2 parents 594af74 + 75caa45
commit 2400d16efc6cfb36554cd64e07c90ba544b2099c
diff --git a/.gitignore b/.gitignore
@@ -1,3 +1,4 @@
 .c9revisions
 *~
 node_modules
+*.pyc
diff --git a/Install/esri2open/__init__.py b/Install/esri2open/__init__.py
@@ -1,2 +1,3 @@
-from esri2open import toOpen, writeFile, closeUp, closeJSON
-from prepare import prepareFile, prepareGeoJSON
+from esri2open import toOpen, writeFile, closeUp, closeJSON,closeTOPOJSON
+from prepare import prepareFile, prepareGeoJSON,prepareTOPO
+from utilities import getExt,getName
diff --git a/Install/esri2open/esri2open.py b/Install/esri2open/esri2open.py
@@ -33,19 +33,27 @@ def closeCSV(out):
     out[1].close()
     return True
 
+def closeTOPOJSON(out):
+    file = open(out['out'],'w')
+    out['topo'].dump(file)
+    file.close()
+    return True
+
 def closeUp(out,fileType):
     if fileType == "geojson":    
         return closeJSON(out)
     elif fileType == "csv":    
         return closeCSV(out)
-    if fileType == "json":    
+    elif fileType == "json":    
         return closeJSON(out)
+    elif fileType == "topojson":    
+        return closeTOPOJSON(out)
     else:
         return False
 
 #this is the meat of the function, we could put it into a seperate file if we wanted
-def writeFile(outArray,featureClass,fileType,includeGeometry, first=True):
-    parser = parse(outArray,featureClass,fileType,includeGeometry,first)
+def writeFile(outArray,featureClass,fileType,includeGeometry, first=True, outName = False):
+    parser = parse(outArray,featureClass,fileType,includeGeometry,first,outName)
     #wrap it in a try so we don't lock the database
     try:
         for row in parser.rows:
@@ -71,7 +79,7 @@ def toOpen(featureClass, outJSON, includeGeometry="geojson"):
         AddMessage("this filetype doesn't make sense")
         return
     #geojson needs geometry
-    if fileType=="geojson":
+    if fileType in ("geojson", "topojson"):
         includeGeometry="geojson"
     elif fileType=="sqlite":
         includeGeometry="well known binary"
diff --git a/Install/esri2open/parseGeometry.py b/Install/esri2open/parseGeometry.py
@@ -206,7 +206,8 @@ def getParseFunc(shpType, geo):
             fun = parseMultiPolygon
         else:
             fun = parseMultiPatch
-    if geo=="geojson":
-        return fun
-    elif geo=="well known text":
-        return getWKTFunc(fun)
+
+    if geo=="well known text":
+        return getWKTFunc(fun)
+    else:
+        return fun
diff --git a/Install/esri2open/parseRow.py b/Install/esri2open/parseRow.py
@@ -1,4 +1,4 @@
-from utilities import listFields, getShp, getOID, statusMessage, parseProp, makeInter
+from utilities import listFields, getShp, getOID, statusMessage, parseProp, makeInter,getName
 from arcpy import SpatialReference, SearchCursor  
 from parseGeometry import getParseFunc
 from json import dump
@@ -7,7 +7,7 @@
 wgs84="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision"
 
 class parse:
-    def __init__(self,outFile,featureClass,fileType,includeGeometry, first=True):
+    def __init__(self,outFile,featureClass,fileType,includeGeometry, first=True, outName=False):
         self.outFile = outFile
         self.fileType = fileType
         #first we set put the local variables we'll need
@@ -33,6 +33,13 @@ def __init__(self,outFile,featureClass,fileType,includeGeometry, first=True):
             self.parse = self.parseJSON
         elif fileType=="sqlite":    
             self.parse = self.parseSqlite
+        elif fileType=="topojson":    
+            self.parse = self.parseTOPOJSON
+            if outName:
+                self.oName=outName
+            else:
+                self.oName = getName(featureClass)
+            self.topo = self.outFile['topo'].object_factory(self.oName)
 
     def cleanUp(self,row):
         del row
@@ -48,6 +55,9 @@ def parseCSV(self,row):
                 fc["geometry"]=self.parseGeo(row.getValue(self.shp))
             except:
                 return
+        for key in fc:
+            if isinstance(fc[key],unicode):
+                fc[key] = fc[key].encode('utf_8')
         self.outFile[0].writerow(fc)
 
     def parseGeoJSON(self,row):
@@ -73,6 +83,23 @@ def parseGeoJSON(self,row):
             self.outFile.write(",")
             dump(fc,self.outFile)
 
+    def parseTOPOJSON(self,row):
+        #more messages
+        self.status.update()
+        fc={"type": "Feature"}
+        if self.parseGeo:
+            try:
+                fc["geometry"]=self.parseGeo(row.getValue(self.shp))
+            except:
+                return
+        else:
+            raise NameError("we need geometry for topojson")
+        fc["id"]=row.getValue(self.oid)
+        fc["properties"]=parseProp(row,self.fields, self.shp)
+        if fc["geometry"]=={}:
+            return
+        self.topo(fc)
+
     def parseJSON(self,row):
         #more messages
         self.status.update()
diff --git a/Install/esri2open/prepare.py b/Install/esri2open/prepare.py
@@ -2,6 +2,7 @@
 from utilities import listFields,getShp, parseFieldType
 from sqlite3 import Connection
 from os.path import splitext, split
+from topojson import Topology
 
 def prepareCSV(outJSON,featureClass,fileType,includeGeometry):
     shp=getShp(featureClass)[0]
@@ -58,6 +59,10 @@ def prepareJSON(outJSON,*args):
     out.write("""{"docs":[""")
     return out
 
+def prepareTOPO(outJSON,*args):
+    topo = Topology()
+    return {"topo":topo,'out':outJSON}
+
 def prepareFile(outJSON,featureClass,fileType,includeGeometry):
     if fileType == "geojson":    
         return prepareGeoJSON(outJSON,featureClass,fileType,includeGeometry)
@@ -67,5 +72,7 @@ def prepareFile(outJSON,featureClass,fileType,includeGeometry):
         return prepareJSON(outJSON,featureClass,fileType,includeGeometry)
     elif fileType == "sqlite": 
         return prepareSqlite(outJSON,featureClass,fileType,includeGeometry)
+    elif fileType == "topojson": 
+        return prepareTOPO(outJSON,featureClass,fileType,includeGeometry)
     else:
         return False
diff --git a/Install/esri2open/topojson/__init__.py b/Install/esri2open/topojson/__init__.py
@@ -0,0 +1 @@
+from topology import Topology
diff --git a/Install/esri2open/topojson/arcs.py b/Install/esri2open/topojson/arcs.py
@@ -0,0 +1,56 @@
+from hashtable import Hashtable
+from os import remove
+import shelve
+from tempfile import mkdtemp
+from hashlib import sha1
+from utils import point_compare
+class Arcs:
+    def __init__(self,Q):
+        self.coincidences = Hashtable(Q * 10)
+        self.arcsByPoint = Hashtable(Q * 10)
+        self.pointsByPoint = Hashtable(Q * 10)
+        self.arc_db_path=mkdtemp()+'/arc_db'
+        self.arcs= shelve.open(self.arc_db_path)
+        #self.arcs={}
+        self.length=0
+        self.storage_path = mkdtemp()+'/db'
+        self.db = shelve.open(self.storage_path)
+        #self.db={}
+    def get_index(self,point):
+        return self.pointsByPoint.get(point)
+    def get_point_arcs(self,point):
+        return self.arcsByPoint.get(point)
+    def coincidence_lines(self,point):
+        return self.coincidences.get(point)
+    def peak(self,point):
+        return self.coincidences.peak(point)
+    def push(self,arc):
+        self.arcs[str(self.length)]=arc
+        self.length+=1
+        return self.length
+    def close(self):
+        #pass
+        self.db.close()
+        remove(self.storage_path)
+        self.arcs.close()
+        remove(self.arc_db_path)
+    def get_hash(self,arc):
+        ourhash = sha1()
+        ourhash.update(str(arc))
+        return ourhash.hexdigest()
+    def check(self,arcs):
+        a0 = arcs[0]
+        a1 = arcs[-1]
+        point = a0 if point_compare(a0, a1) < 0 else a1
+        point_arcs = self.get_point_arcs(point)
+        h = self.get_hash(arcs)
+        if h in self.db:
+            return int(self.db[h])
+        else:
+            index = self.length
+            point_arcs.append(arcs)
+            self.db[h]=index
+            self.db[self.get_hash(list(reversed(arcs)))]=~index
+            self.push(arcs)
+            return index
+        
diff --git a/Install/esri2open/topojson/bounds.py b/Install/esri2open/topojson/bounds.py
@@ -0,0 +1,20 @@
+from mytypes import Types
+class Bounds(Types):
+    def __init__(self):
+        self.x0=self.y0=float('inf')
+        self.x1=self.y1=-float('inf')
+    def point (self,point):
+        x = point[0]
+        y = point[1]
+        if x < self.x0:
+            self.x0 = x
+        if x > self.x1:
+            self.x1 = x
+        if y < self.y0:
+            self.y0 = y
+        if y > self.y1:
+            self.y1 = y
+def bound(objects):
+    b=Bounds()
+    b.obj(objects)
+    return b
diff --git a/Install/esri2open/topojson/clockwise.py b/Install/esri2open/topojson/clockwise.py
@@ -0,0 +1,25 @@
+class Clock:
+    def __init__(self,area):
+        self.area=area
+    def clock(self,feature):
+        if 'geometries' in feature:
+            feature['geometries'] = map(self.clock_geometry,feature['geometries'])
+        elif 'geometry' in feature:
+            feature['geometry']=self.clock_geometry(feature['geometry'])
+        return feature
+    def clock_geometry(self,geo):
+        if 'type' in geo:
+            if geo['type']=='Polygon' or geo['type']=='MultiLineString':
+                geo['coordinates'] = self.clockwise_polygon(geo['coordinates'])
+            elif geo['type']=='MultiPolygon':
+                geo['coordinates'] = map(lambda x:self.clockwise_polygon(x),geo['coordinates'])
+            elif geo['type']=='LineString':
+                geo['coordinates'] = self.clockwise_ring(geo['coordinates'])
+        return geo
+    def clockwise_polygon(self,rings):
+        return map(lambda x:self.clockwise_ring(x),rings)
+    def clockwise_ring(self,ring):
+        if self.area(ring) > 0:
+            return list(reversed(ring))
+        else:
+            return ring
diff --git a/Install/esri2open/topojson/coordinatesystems.py b/Install/esri2open/topojson/coordinatesystems.py
@@ -0,0 +1,118 @@
+# coding=utf8
+from math import sqrt, pi, cos, sin, atan2, tan, atan, asin
+
+PI4 = pi / 4
+RADIANS = pi / 180
+
+
+class BaseCoordinateSystem(object):
+    name = "BaseCoordinateSystem"
+
+    def format_distance(self, distance):
+        return str(distance)
+
+    def ring_area(self, ring):
+        raise Exception('Not implemented')
+
+    def triangle_area(self, triangle):
+        raise Exception('Not implemented')
+
+    def distance(self, x0, y0, x1, y1):
+        raise Exception('Not implemented')
+
+    def absolute_area(self, area):
+        return abs(area)
+
+
+class Cartesian(BaseCoordinateSystem):
+    name = "cartesian"
+
+    def ring_area(self, ring):
+        calc_area = lambda p1, p2: p1[1] * p2[0] - p1[0] * p2[1]
+        # last and first
+        area = calc_area(ring[-1], ring[0])
+        for i, p in enumerate(ring[1:]):   # skip first so p is current and
+            area += calc_area(p, ring[i])  # ring[i] is the previous
+        return area * 0.5
+
+    def triangle_area(self, triangle):
+        return abs(
+            (triangle[0][0] - triangle[2][0]) * (triangle[1][1] - triangle[0][1]) -
+            (triangle[0][0] - triangle[1][0]) * (triangle[2][1] - triangle[0][1])
+        )
+
+    def distance(self, x0, y0, x1, y1):
+        dx = x0 - x1
+        dy = y0 - y1
+        return sqrt(dx * dx + dy * dy)
+
+
+class Spherical(BaseCoordinateSystem):
+    name = 'spherical'
+
+    def haversin(self, x):
+        return sin(x / 2) ** 2
+
+    def format_distance(self, distance):
+        km = distance * 6371.0
+        if km > 1:
+            return u"{:0.03f}km".format(km)
+        else:
+            return u"{:0.03f} ({0.03f}°)".format(km * 1000, distance * 180 / pi)
+
+    def ring_area(self, ring):
+        if len(ring) == 0:
+            return 0
+        area = 0
+        p = ring[0]
+        PI4 = 1
+        lambda_ = p[0] * RADIANS
+        phi = p[1] * RADIANS / 2.0 + PI4
+        lambda0 = lambda_
+        cosphi0 = cos(phi)
+        sinphi0 = sin(phi)
+        for pp in ring[1:]:
+            lambda_ = pp[0] * RADIANS
+            phi = pp[1] * RADIANS / 2.0 + PI4
+            # Spherical excess E for a spherical triangle with vertices: south pole,
+            # previous point, current point.  Uses a formula derived from Cagnoli’s
+            # theorem.  See Todhunter, Spherical Trig. (1871), Sec. 103, Eq. (2).
+            dlambda = lambda_ - lambda0
+            cosphi = cos(phi)
+            sinphi = sin(phi)
+            k = sinphi0 * sinphi
+            u = cosphi0 * cosphi + k * cos(dlambda)
+            v = k * sin(dlambda)
+            area += atan2(v, u)
+            #Advance the previous point.
+            lambda0 = lambda_
+            cosphi0 = cosphi
+            sinphi0 = sinphi
+        return 2 * area
+
+    def absolute_area(self, area):
+        return area + 4 * pi if area < 0 else area
+
+    def triangle_area(self, triangle):
+        def distance(a, b):    # why 2 implementations? I don't know, original has the same question in comments
+            x0, y0, x1, y1 = [(n * RADIANS) for n in (a + b)]
+            delta_lambda = x1 - x0
+            return atan2(
+                sqrt(
+                    (cos(x1) * sin(delta_lambda)) ** 2 +
+                    (cos(x0) * sin(x1) - sin(x0) * cos(x1) * cos(delta_lambda)) ** 2
+                ),
+                sin(x0) * sin(x1) + cos(x0) * cos(x1) * cos(delta_lambda)
+            )
+        a = distance(triangle[0], triangle[1])
+        b = distance(triangle[1], triangle[2])
+        c = distance(triangle[2], triangle[0])
+        s = (a + b + c) / 2.0
+        return 4 * atan(sqrt(max(0, tan(s / 2.0) * tan((s - a) / 2.0) * tan((s - b) / 2.0) * tan((s - c) / 2.0))))
+
+    def distance(self, x0, y0, x1, y1):
+        x0, y0, x1, y1 = [(n * RADIANS) for n in [x0, y0, x1, y1]]
+        return 2.0 * asin(sqrt(self.haversin(y1 - y0) + cos(y0) * cos(y1) * self.haversin(x1 - x0)))
+
+
+systems = {'cartesian': Cartesian(), 'spherical': Spherical()}