Wkt parsing into tuples and args now works for easier parsing. Next just need to make them into objs.

pycrs/datums.py

@@ -23,6 +23,13 @@ class NAD83:
     ellipsdef = "" # ellipsoids.WGS72()
     to_wgs84 = None
 
+class NAD27:
+    proj4 = "NAD27"
+    ogc_wkt = "D_North_American_1927"
+    esri_wkt = "D_North_American_1927"
+    
+    ellipsdef = "" # ellipsoids...
+    to_wgs84 = None
 
 class Unknown:
     proj4 = "unknown" # no datum name, just ellips + towgs84 params...

pycrs/ellipsoids.py

@@ -32,3 +32,22 @@ class GRS80:
 
     semimaj_ax = 6378137.0
     inv_flat = 298.257222101
+
+class Clarke1866:
+    proj4 = "clrk66"
+    ogc_wkt = "Clarke_1866"
+    esri_wkt = "Clarke_1866"
+
+    semimaj_ax = 6378206.4
+    inv_flat = 294.9786982
+
+class Airy1830:
+    proj4 = "airy"
+    ogc_wkt = "Airy 1830"
+    esri_wkt = "Airy_1830"
+
+    semimaj_ax = 6377563.396
+    inv_flat = 299.3249646
+
+
+    

pycrs/parameters.py

@@ -444,7 +444,6 @@ class FalseEasting:
         return "+x_0=%s" %self.value
 
     def to_ogc_wkt(self):
-        # the stuff that comes after UNITS[... # must be combined with metermultiplier in a unit class to make wkt
         return 'PARAMETER["false_easting", %s]' % self.value
 
     def to_esri_wkt(self):
@@ -464,11 +463,39 @@ class FalseNorthing:
         return "+y_0=%s" %self.value
 
     def to_ogc_wkt(self):
-        # the stuff that comes after UNITS[... # must be combined with metermultiplier in a unit class to make wkt
         return 'PARAMETER["false_northing", %s]' % self.value
 
     def to_esri_wkt(self):
         return self.to_ogc_wkt()
+
+##+h       Satellite height
+class SatelliteHeight:
+    def __init__(self, value):
+        self.value = value
+
+    def to_proj4(self):
+        return "+h=%s" %self.value
+
+    def to_ogc_wkt(self):
+        return 'PARAMETER["satellite_height", %s]' % self.value
+
+    def to_esri_wkt(self):
+        return self.to_ogc_wkt()
+
+##+tilt     Tilt angle
+class TiltAngle:
+    def __init__(self, value):
+        self.value = value
+
+    def to_proj4(self):
+        return "+tilt=%s" %self.value
+
+    def to_ogc_wkt(self):
+        raise Exception("Paramater not supported by OGC WKT")
+
+    def to_esri_wkt(self):
+        raise Exception("Paramater not supported by ESRI WKT")
+
     
 # then the final CRS object which is instantiated with all of these?
 # remember to use +no_defs when outputting to proj4

pycrs/parser.py

@@ -3,6 +3,7 @@
 # possible use module: https://github.com/rockdoc/grabbag/wiki/CRS-WKT-Parser
 # also note some paramter descriptions: http://www.geoapi.org/3.0/javadoc/org/opengis/referencing/doc-files/WKT.html
 # and see gdal source code: http://gis.stackexchange.com/questions/129764/how-are-esri-wkt-projections-different-from-ogc-wkt-projections
+# especially: http://fossies.org/windows/misc/saga_2.1.4_x64.zip/saga_2.1.4_x64/saga_prj.dic
 
 from . import datums
 from . import ellipsoids
@@ -36,12 +37,133 @@ def from_sr_code(code):
 ##    # parse arguments into components
 ##    # use args to create crs
 ##    pass
-##        
-##def from_ogc_wkt(string):
-##    # parse arguments into components
-##    # use args to create crs
-##    pass
-##
+        
+def from_ogc_wkt(string):
+    # remove newlines and multi spaces
+    string = " ".join(string.split())
+    
+    # parse arguments into components
+    def _consume_bracket(chars, char):
+        "char must be the opening bracket"
+        consumed = ""
+        depth = 1
+        while char and depth > 0:
+            consumed += char
+            char = next(chars, None)
+            # update depth level
+            if char == "[":
+                depth += 1
+            elif char == "]":
+                depth -= 1
+        consumed += char # consume the last closing char too
+        return consumed
+    
+    def _consume_quote(chars, char, quotechar):
+        "char and quotechar must be the opening quote char"
+        consumed = ""
+        # consume the first opening char
+        consumed += char 
+        char = next(chars, None)
+        # consume inside
+        while char and char != quotechar:
+            consumed += char
+            char = next(chars, None)
+        # consume the last closing char too
+        consumed += char 
+        return consumed
+    
+    def _next_elem(chars, char):
+        "char can be any char"
+        header = ""
+        # skip until next header
+        while not char.isalpha():
+            char = next(chars, None)
+        # first consume the element text header
+        while char.isalpha():
+            header += char
+            char = next(chars, None)
+        # skip until next brackets (in case of spaces)
+        while char != "[":
+            char = next(chars, None)
+        # then consume the element bracket contents
+        if char == "[":
+            content = _consume_bracket(chars, char)
+        char = next(chars, None)
+        # split content into args list
+        content = content[1:-1] # remove enclosing brackets
+        content = _split_except(content)
+        # recursively load all subelems
+        for i,item in enumerate(content):
+            if isinstance(item, str) and "[" in item:
+                chars = (char for char in item)
+                char = next(chars)
+                item = _next_elem(chars, char)
+                content[i] = item
+        return header, content
+    
+    def _clean_value(string):
+        string = string.strip()
+        try: string = float(string)
+        except: pass
+        return string
+    
+    def _split_except(string):
+        chars = (char for char in string)
+        char = next(chars)
+        items = []
+        consumed = ""
+        while char:
+            # dont split on quotes, just consume it
+            if char in ("'", '"'):
+                consumed += _consume_quote(chars, char, char)
+            # dont split inside brackets, just consume it
+            elif char == "[":
+                consumed += _consume_bracket(chars, char)
+            # make new splititem
+            elif char == ",":
+                consumed = _clean_value(consumed)
+                items.append(consumed)
+                consumed = ""
+            # consume normal char
+            elif char:
+                consumed += char
+            # next
+            char = next(chars, None)
+        # append last item too
+        consumed = _clean_value(consumed)
+        items.append(consumed)
+        return items
+
+    # load into nested tuples and arglists
+    crstuples = []
+    chars = (char for char in string)
+    char = next(chars)
+    while char:
+        header,content = _next_elem(chars, char)
+        crstuples.append((header, content))
+        char = next(chars, None)
+
+    # parse into actual crs objects
+    for header, content in crstuples:
+        # toplevel collection
+        if header.upper() == "PROJCS":
+            # find name
+            # find geogcs elem
+            # find projection elem
+            # find params
+            # find unit
+            pass
+            #projcs = parameters.ProjCS("Unknown", geogcs, proj, params, unit)
+
+        elif header.upper() == "GEOGCS":
+            pass
+
+            
+
+        
+    # use args to create crs
+    return crstuples
+
 ##def from_unknown_wkt(string):
 ##    # detect if ogc wkt or esri wkt
 ##    # TIPS: esri wkt datums all use "D_" before the datum name
@@ -61,6 +183,13 @@ def from_proj4(string):
 
     # INIT CODES...?
     # eg, +init=epsg code...?
+    if "+init" in partdict:
+
+        codetype, code = partdict["+init"].split(":")
+        if codetype == "EPSG":
+            crs = from_epsg_code(code)
+            # need to get individual elements from crs
+            # maybe from web...
 
     # DATUM
 
@@ -68,11 +197,15 @@ def from_proj4(string):
     if "+datum" in partdict:
         
         # get predefined datum def
-        if partdict["+datum"] == "WGS84":
-            datumdef = datums.WGS84()
-
-        elif partdict["+datum"] == "NAD83":
-            datumdef = datums.NAD83()
+        for itemname in dir(datums):
+            item = getattr(datums, itemname)
+            try:
+                item = item()
+                if hasattr(item, "proj4") and partdict["+datum"] == item.proj4:
+                    datumdef = item
+                    break
+            except:
+                pass
 
         else:
             datumdef = datums.Unknown()
@@ -86,17 +219,15 @@ def from_proj4(string):
     if "+ellps" in partdict:
 
         # get predefined ellips def
-        if partdict["+ellps"] == "WGS84":
-            ellipsdef = ellipsoids.WGS84()
-
-        elif partdict["+ellps"] == "WGS72":
-            ellipsdef = ellipsoids.WGS72()
-
-        elif partdict["+ellps"] == "intl":
-            ellipsdef = ellipsoids.International()
-
-        elif partdict["+ellps"] == "GRS80":
-            ellipsdef = ellipsoids.GRS80()
+        for itemname in dir(ellipsoids):
+            item = getattr(ellipsoids, itemname)
+            try:
+                item = item()
+                if hasattr(item, "proj4") and partdict["+ellps"] == item.proj4:
+                    ellipsdef = item
+                    break
+            except:
+                pass
 
     else:
         raise Exception("Could not find required +ellps element")
@@ -163,22 +294,15 @@ def from_proj4(string):
     if "+proj" in partdict:
 
         # get predefined proj def
-        if partdict["+proj"] == "robin":
-            projdef = projections.Robinson()
-
-        elif partdict["+proj"] == "omerc":
-            projdef = projections.ObliqueMercator()
-
-        elif partdict["+proj"] == "longlat":
-            projdef = None
-            # set geogcs axis in correct order
-
-        elif partdict["+proj"] == "latlong":
-            projdef = None
-            # set geogcs axis in correct order
-
-            # ALSO SHOULDNT EXCLUDE +proj, NEED WAY TO INCLUDE IT...
-            # ...
+        for itemname in dir(projections):
+            item = getattr(projections, itemname)
+            try:
+                item = item()
+                if hasattr(item, "proj4") and partdict["+proj"] == item.proj4:
+                    projdef = item
+                    break
+            except:
+                pass
 
         else:
             projdef = None
@@ -278,6 +402,12 @@ def from_proj4(string):
         ## create unitobj
         unit = parameters.Unit(unittype, metmulti)
 
+        # SATELLITE HEIGHT
+        if "+h" in partdict:
+            val = partdict["+h"]
+            obj = parameters.SatelliteHeight(val)
+            params.append(obj)
+
         # PROJCS
 
         projcs = parameters.ProjCS("Unknown", geogcs, proj, params, unit)

pycrs/projections.py

@@ -12,3 +12,114 @@ class ObliqueMercator:
     proj4 = "omerc"
     ogc_wkt = "Hotine_Oblique_Mercator"
     
+class AlbersEqualArea:
+    proj4 = "aea"
+    ogc_wkt = "Albers"
+
+class CylindricalEqualArea:
+    proj4 = "cea"
+    ogc_wkt = "Cylindrical_Equal_Area"
+    esri_wkt = "Cylindrical_Equal_Area"
+    
+class EquiDistantConic:
+    proj4 = "eqdc"
+    ogc_wkt = "Equidistant_Conic"
+
+class TransverseMercator:
+    proj4 = "tmerc"
+    ogc_wkt = "Transverse_Mercator"
+
+class GallStereographic:
+    proj4 = "gall"
+    ogc_wkt = "Gall_Stereographic"
+
+class Gnomonic:
+    proj4 = "gnom"
+    ogc_wkt = "Gnomonic"
+
+class LambertAzimuthalEqualArea:
+    proj4 = "laea"
+    ogc_wkt = "Lambert_Azimuthal_Equal_Area"
+
+class MillerCylindrical:
+    proj4 = "mill"
+    ogc_wkt = "Miller_Cylindrical"
+
+class Mollweide:
+    proj4 = "moll"
+    ogc_wkt = "Mollweide"
+
+class ObliqueStereographic:
+    proj4 = "sterea"
+    ogc_wkt = "Oblique_Stereographic"
+
+class Orthographic:
+    proj4 = "ortho"
+    ogc_wkt = "Orthographic"
+
+class PolarStereographic:
+    proj4 = "stere"
+    ogc_wkt = "Polar_Stereographic"
+
+class Sinusoidal:
+    proj4 = "sinu"
+    ogc_wkt = "Sinusoidal"
+
+class VanDerGrinten:
+    proj4 = "vandg"
+    ogc_wkt = "VanDerGrinten"
+
+class Equirectangular:
+    proj4 = "eqc"
+    ogc_wkt = "Equirectangular"
+
+class LambertConformalConic:
+    proj4 = "lcc"
+    ogc_wkt = "Lambert_Conformal_Conic"
+
+class Krovak:
+    proj4 = "krovak"
+    ogc_wkt = "Krovak"
+
+class NearSidedPerspective:
+    proj4 = "nsper"
+    ogc_wkt = "Near_sided_perspective"
+
+class TiltedPerspective:
+    proj4 = "tsper"
+    ogc_wkt = "Tilted_perspective"
+
+class InteruptedGoodeHomolosine:
+    proj4 = "igh"
+    ogc_wkt = "Interrupted_Goodes_Homolosine"
+
+class Larrivee:
+    proj4 = "larr"
+    ogc_wkt = "Larrivee"
+
+class LamberEqualAreaConic:
+    proj4 = "leac"
+    ogc_wkt = "Lambert_Equal_Area_Conic"
+
+class Mercator:
+    proj4 = "merc"
+    ogc_wkt = "Mercator"
+
+class ObliqueCylindricalEqualArea:
+    proj4 = "ocea"
+    ogc_wkt = "Oblique_Cylindrical_Equal_Area"
+
+        
+
+
+
+
+
+
+
+
+    
+
+    
+
+

pycrs/units.py

@@ -9,3 +9,8 @@ class Degree:
     proj4 = "degrees"
     ogc_wkt = "degree"
     esri_wkt = "Degree"
+
+class Feet:
+    proj4 = "..."
+    ogc_wkt = "Foot_US"
+    esri_wkt = "Foot_US"

tester.py

@@ -1,8 +1,24 @@
 import pycrs
 
+###########################
+# OGC WKT
+
+crs = pycrs.parser.from_esri_code(54030)
+wkt = crs.to_ogc_wkt()
+print pycrs.parser.from_ogc_wkt(wkt)
+
+
+
+dssdfsd
+
+###########################
+# PROJ4
+
 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"
 
 #crs = pycrs.parser.from_esri_code(54030)
 #crs = pycrs.parser.from_sr_code(7898)
@@ -21,6 +37,7 @@ 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()
@@ -35,16 +52,22 @@ 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))
+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:
-    if feat.geometry.type == "Polygon":
-        try: c.draw_geojson(feat.geometry, outlinecolor="white")
-        except:
-            # NOTE: feat.__geo_interface__ is one level too high maybe??
-            print feat.geometry
+    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()