Basic ogc wkt parsing into crsobj works.

Next:
- Comb through to smooth out kinks, like how to handle default values
such as datum (not always needed) or unit when not specified (on parsing
or on outputting). Or handling basic lat long crs without proj. Adding
wkt axis parsing. Cleaning messy wkt parsing maybe.
- And finally add esri wkt, copy paste mostly, but some small
differences.

pycrs/loader.py

@@ -33,18 +33,25 @@ def from_file(filepath):
         return parser.from_esri_wkt(string)
     
     elif filepath.endswith((".geojson",".json")):
-        crsinfo = json.load(filepath)["crs"]
+        raw = open(filepath).read()
-        
+        geoj = json.loads(raw)
-        if crsinfo["type"] == "name":
+        if "crs" in geoj:
-            string = crsinfo["properties"]["name"]
+            crsinfo = geoj["crs"]
-            return parser.from_unknown_text(string)
             
-        elif crsinfo["type"] == "link":
+            if crsinfo["type"] == "name":
-            url = crsinfo["properties"]["name"]
+                string = crsinfo["properties"]["name"]
-            type = crsinfo["properties"].get("type")
+                return parser.from_unknown_text(string)
-            return from_url(url, format=type)
+                
-            
+            elif crsinfo["type"] == "link":
-        else: raise Exception("invalid geojson crs type: must be either name or link")
+                url = crsinfo["properties"]["name"]
+                type = crsinfo["properties"].get("type")
+                return loader.from_url(url, format=type)
+                
+            else: raise Exception("invalid geojson crs type: must be either name or link")
+
+        else:
+            # assume default wgs84 as per the spec
+            return parser.from_epsg_code("4326")
 
 ##    elif filepath.endswith((".tif",".tiff",".geotiff")):
 ##        pass

pycrs/parameters.py

@@ -163,6 +163,7 @@ class ProjCS:
         self.geogcs = geogcs
         self.proj = proj
         self.params = params
+        self.unit = unit
         if twin_ax == None:
             # default axes
             twin_ax = directions.East(), directions.North()
@@ -171,6 +172,7 @@ class ProjCS:
     def to_proj4(self):
         string = "%s %s " % (self.proj.to_proj4(), self.geogcs.to_proj4())
         string += " ".join(param.to_proj4() for param in self.params)
+        string += " %s" % self.unit.to_proj4()
         # in proj4, axis only applies to the cs, ie the projcs (not the geogcs, where wkt can specify with axis)
         string += " +axis=" + self.twin_ax[0].proj4 + self.twin_ax[1].proj4 + "u" # up set as default because only proj4 can set it I think...
         return string
@@ -178,6 +180,7 @@ class ProjCS:
     def to_ogc_wkt(self):
         string = 'PROJCS["%s", %s, %s, ' % (self.name, self.geogcs.to_ogc_wkt(), self.proj.to_ogc_wkt() )
         string += ", ".join(param.to_ogc_wkt() for param in self.params)
+        string += ', %s' % self.unit.to_ogc_wkt()
         string += ', AXIS["X", %s], AXIS["Y", %s]]' % (self.twin_ax[0].ogc_wkt, self.twin_ax[1].ogc_wkt )
         return string
     
@@ -356,10 +359,10 @@ class DatumShift:
         self.value = value
 
     def to_proj4(self):
-        return "+towgs84=%s" %",".join((str(val) for val in self.value))
+        return "+towgs84=%s" %",".join((bytes(val) for val in self.value))
 
     def to_ogc_wkt(self):
-        return "TOWGS84[%s]" %",".join((str(val) for val in self.value))
+        return "TOWGS84[%s]" %",".join((bytes(val) for val in self.value))
 
     def to_esri_wkt(self):
         raise Exception("Paramater not supported by ESRI WKT")
@@ -374,7 +377,7 @@ class MeterMultiplier:
 
     def to_ogc_wkt(self):
         # the stuff that comes after UNITS["meter", ... # must be combined with unittype in a unit class to make wkt
-        return str(self.value)
+        return bytes(self.value)
 
     def to_esri_wkt(self):
         return self.to_ogc_wkt()
@@ -394,7 +397,7 @@ class UnitType:
 
     def to_ogc_wkt(self):
         # the stuff that comes after UNITS[... # must be combined with metermultiplier in a unit class to make wkt
-        return str(self.value.ogc_wkt)
+        return bytes(self.value.ogc_wkt)
 
     def to_esri_wkt(self):
         return self.to_ogc_wkt()

pycrs/parser.py

@@ -73,7 +73,7 @@ def from_ogc_wkt(string):
         return consumed
     
     def _next_elem(chars, char):
-        "char can be any char"
+        "char must be the first char of the text that precedes brackets"
         header = ""
         # skip until next header
         while not char.isalpha():
@@ -144,25 +144,150 @@ def from_ogc_wkt(string):
         char = next(chars, None)
 
     # parse into actual crs objects
-    for header, content in crstuples:
+    def _parse(header, content):
-        # toplevel collection
         if header.upper() == "PROJCS":
             # find name
-            # find geogcs elem
+            name = content[0].strip('"')
+            # find geogcs elem (by running parse again)
+            subheader, subcontent = content[1]
+            geogcs = _parse(subheader, subcontent)
             # find projection elem
+            subheader, subcontent = content[2]
+            for itemname in dir(projections):
+                item = getattr(projections, itemname)
+                try:
+                    item = item()
+                    if hasattr(item, "ogc_wkt") and subcontent[0].strip('"') == item.ogc_wkt:
+                        projdef = item
+                        break
+                except:
+                    pass
+            else:
+                projdef = None
+            proj = parameters.Projection(projdef)
             # find params
+            params = []
+            for part in content:
+                if isinstance(part, tuple):
+                    subheader,subcontent = part
+                    if subheader == "PARAMETER":
+                        if subcontent[0] == '"Azimuth"':
+                            item = parameters.Azimuth(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"scale_factor"':
+                            item = parameters.ScalingFactor(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"latitude_of_origin"':
+                            item = parameters.LatitudeOrigin(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"standard_parallel_1"':
+                            item = parameters.LatitudeFirstStndParallel(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"standard_parallel_2"':
+                            item = parameters.LatitudeSecondStndParallel(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"Standard_Parallel_1"':
+                            item = parameters.LatitudeTrueScale(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"Central_Meridian"':
+                            item = parameters.CentralMeridian(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"Longitude_Of_Center"':
+                            item = parameters.LongitudeCenter(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"false_easting"':
+                            item = parameters.FalseEasting(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"false_northing"':
+                            item = parameters.FalseNorthing(subcontent[1])
+                            params.append(item)
+                        elif subcontent[0] == '"satellite_height"':
+                            item = parameters.SatelliteHeight(subcontent[1])
+                            params.append(item)
             # find unit
-            pass
+            for part in content:
-            #projcs = parameters.ProjCS("Unknown", geogcs, proj, params, unit)
+                if isinstance(part, tuple):
+                    subheader,subcontent = part
+                    if subheader == "UNIT":
+                        break
+            if subcontent[0].strip('"') == "Meters":
+                unittype = parameters.UnitType(units.Meter())
+            elif subcontent[0].strip('"') == "degree":
+                unittype = parameters.UnitType(units.Degree())
+            metmult = parameters.MeterMultiplier(subcontent[1])
+            unit = parameters.Unit(unittype, metmult)
+            # find twin axis maybe
+##            if len(content) >= 6:
+##                twinax = (parameters.Axis(
+##            else:
+##                twinax = None
+            # put it all together
+            projcs = parameters.ProjCS("Unknown", geogcs, proj, params, unit) #, twinax)
+            return projcs
 
         elif header.upper() == "GEOGCS":
-            pass
+            # name
-
+            name = content[0]
-            
+            # datum
+            subheader, subcontent = content[1]
+            ## datum name
+            for itemname in dir(datums):
+                item = getattr(datums, itemname)
+                try:
+                    item = item()
+                    if hasattr(item, "ogc_wkt") and subcontent[0].strip('"') == item.ogc_wkt:
+                        datumdef = item
+                        break
+                except:
+                    pass
+            else:
+                datumdef = datums.Unknown()
+            ## datum ellipsoid
+            subsubheader, subsubcontent = subcontent[1]
+            for itemname in dir(ellipsoids):
+                item = getattr(ellipsoids, itemname)
+                try:
+                    item = item()
+                    if hasattr(item, "ogc_wkt") and subsubcontent[0].strip('"') == item.ogc_wkt:
+                        ellipsdef = item
+                        break
+                except:
+                    pass
+            else:
+                ellipsdef = None
+            ellipsoid = parameters.Ellipsoid(ellipsdef, subsubcontent[1], subsubcontent[2])
+            ## datum shift
+            if len(subcontent) >= 3:
+                subsubheader, subsubcontent = subcontent[2]
+                datumshift = parameters.DatumShift(subsubcontent)
+            else:
+                datumshift = None
+            ## put it all togehter
+            datum = parameters.Datum(datumdef, ellipsoid, datumshift)
+            # prime mer
+            subheader, subcontent = content[2]
+            prime_mer = parameters.PrimeMeridian(subcontent[1])
+            # angunit
+            subheader, subcontent = content[3]
+            if subcontent[0].strip('"') == "Meters":
+                unittype = parameters.UnitType(units.Meter())
+            elif subcontent[0].strip('"') == "degree":
+                unittype = parameters.UnitType(units.Degree())
+            metmult = parameters.MeterMultiplier(subcontent[1])
+            angunit = parameters.AngularUnit(unittype, metmult)
+            # twin axis
+            # ...
+            # put it all together
+            geogcs = parameters.GeogCS(name, datum, prime_mer, angunit, twin_ax=None)
+            return geogcs
 
+    # toplevel collection
+    header, content in crstuples[0]
+    toplevel = _parse(header, content)
+    crs = parameters.CRS(toplevel)
         
     # use args to create crs
-    return crstuples
+    return crs
 
 ##def from_unknown_wkt(string):
 ##    # detect if ogc wkt or esri wkt
@@ -408,6 +533,12 @@ def from_proj4(string):
             obj = parameters.SatelliteHeight(val)
             params.append(obj)
 
+        # TILT ANGLE
+        if "+tilt" in partdict:
+            val = partdict["+tilt"]
+            obj = parameters.TiltAngle(val)
+            params.append(obj)
+
         # PROJCS
 
         projcs = parameters.ProjCS("Unknown", geogcs, proj, params, unit)

tester.py

@@ -1,73 +1,99 @@
 import pycrs
 
+
+
+
+###########################
+# Drawing routine for testing
+def render_world(crs):
+    import urllib2
+    import json
+    import pygeoj
+    import pyagg
+    import pyproj
+    import random
+    c = pyagg.Canvas(600,600)
+    c.geographic_space()
+    raw = urllib2.urlopen("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json").read()
+    rawdict = json.loads(raw)
+    data = pygeoj.load(data=rawdict)
+    proj4 = crs.to_proj4()
+    fromproj = pyproj.Proj("+init=EPSG:4326")
+    toproj = pyproj.Proj(proj4)
+
+    print data.bbox
+    for feat in data:
+        if feat.geometry.type == "Polygon":
+            feat.geometry.coordinates = [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1]))
+                                         for ring in feat.geometry.coordinates]
+        elif feat.geometry.type == "MultiPolygon":
+            feat.geometry.coordinates = [
+                                        [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1]))
+                                         for ring in poly]
+                                         for poly in feat.geometry.coordinates]
+    data.add_all_bboxes()
+    xmins, ymins, xmaxs, ymaxs = zip(*(feat.geometry.bbox for feat in data))
+    bbox = (min(xmins), min(ymins), max(xmaxs), max(ymaxs))#33333, 33333) # to avoid inf bounds and no render
+    # print data.bbox #PYGEOJ data.bbox doesnt update.........FIX!
+    print bbox
+    c.zoom_bbox(*bbox)
+    for feat in data:
+        try: c.draw_geojson(feat.geometry,
+                            fillcolor=tuple(random.randrange(255) for _ in range(3)),
+                            outlinecolor="white")
+        except:
+            # NOTE: feat.__geo_interface__ is one level too high maybe??
+            print feat.geometry
+    c.view()
+
+
+
+
+
+
 ###########################
 # OGC WKT
-
+print "testing ogc wkt"
-crs = pycrs.parser.from_esri_code(54030)
+crs = pycrs.parser.from_sr_code(54030)
+#crs = pycrs.parser.from_sr_code(22)
+print crs.to_proj4()
 wkt = crs.to_ogc_wkt()
-print pycrs.parser.from_ogc_wkt(wkt)
+print "Original:", wkt
+print "Reconstructed:", pycrs.parser.from_ogc_wkt(wkt).to_ogc_wkt()
+
+render_world(crs)
 
 
 
-dssdfsd
 
 ###########################
 # PROJ4
-
+print "testing proj4"
-proj4 = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
+#proj4 = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
 #proj4 = "+proj=longlat +ellps=WGS84 +datum=WGS84"
 #proj4 = "+proj=aea +lat_1=24 +lat_2=31.5 +lat_0=24 +lon_0=-84 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs "
 proj4 = "+proj=larr +datum=WGS84 +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
-proj4 = "+proj=nsper +datum=WGS84 +ellps=WGS84 +lon_0=44 +lat_0=30 +h=200000000"
+#proj4 = "+proj=nsper +datum=WGS84 +ellps=WGS84 +lon_0=44 +lat_0=30 +h=200000000"
 
 #crs = pycrs.parser.from_esri_code(54030)
 #crs = pycrs.parser.from_sr_code(7898)
 #crs = pycrs.parser.from_epsg_code(4324)
 #crs = pycrs.parser.from_sr_code(6618)
 crs = pycrs.parser.from_proj4(proj4)
+#crs = pycrs.parser.from_ogc_wkt(wkt)
 
 print crs.to_proj4()
 
+render_world(crs)
 
 
 
 
-import urllib2
+###########################
-import json
+# PRJ FILE
-import pygeoj
+crs = pycrs.loader.from_file("testfiles/natearth.prj")
-import pyagg
+print crs
-import pyproj
+
-import random
+
-c = pyagg.Canvas(600,600)
+
-c.geographic_space()
-raw = urllib2.urlopen("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json").read()
-rawdict = json.loads(raw)
-data = pygeoj.load(data=rawdict)
-proj4 = crs.to_proj4()
-fromproj = pyproj.Proj("+init=EPSG:4326")
-toproj = pyproj.Proj(proj4)
 
-print data.bbox
-for feat in data:
-    if feat.geometry.type == "Polygon":
-        feat.geometry.coordinates = [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1]))
-                                     for ring in feat.geometry.coordinates]
-    elif feat.geometry.type == "MultiPolygon":
-        feat.geometry.coordinates = [
-                                    [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1]))
-                                     for ring in poly]
-                                     for poly in feat.geometry.coordinates]
-data.add_all_bboxes()
-xmins, ymins, xmaxs, ymaxs = zip(*(feat.geometry.bbox for feat in data))
-bbox = (min(xmins), min(ymins), 33333, 33333) # to avoid inf bounds and no render
-# print data.bbox #PYGEOJ data.bbox doesnt update.........FIX!
-print bbox
-c.zoom_bbox(*bbox)
-for feat in data:
-    try: c.draw_geojson(feat.geometry,
-                        fillcolor=tuple(random.randrange(255) for _ in range(3)),
-                        outlinecolor="white")
-    except:
-        # NOTE: feat.__geo_interface__ is one level too high maybe??
-        print feat.geometry
-c.view()

testfiles/natearth.prj

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]