diff --git a/pycrs/datums.py b/pycrs/datums.py index ac77fee..20a1606 100644 --- a/pycrs/datums.py +++ b/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... diff --git a/pycrs/datums.pyc b/pycrs/datums.pyc index bdf7a95..cdade30 100644 Binary files a/pycrs/datums.pyc and b/pycrs/datums.pyc differ diff --git a/pycrs/ellipsoids.py b/pycrs/ellipsoids.py index 101bf89..37ad624 100644 --- a/pycrs/ellipsoids.py +++ b/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 + + + diff --git a/pycrs/ellipsoids.pyc b/pycrs/ellipsoids.pyc index 0249cb0..073532d 100644 Binary files a/pycrs/ellipsoids.pyc and b/pycrs/ellipsoids.pyc differ diff --git a/pycrs/parameters.py b/pycrs/parameters.py index 6bdd155..cfe06e6 100644 --- a/pycrs/parameters.py +++ b/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 diff --git a/pycrs/parameters.pyc b/pycrs/parameters.pyc index 952d451..f678fd6 100644 Binary files a/pycrs/parameters.pyc and b/pycrs/parameters.pyc differ diff --git a/pycrs/parser.py b/pycrs/parser.py index 311da53..f1fa47f 100644 --- a/pycrs/parser.py +++ b/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) diff --git a/pycrs/parser.pyc b/pycrs/parser.pyc index ca4af78..4a3a414 100644 Binary files a/pycrs/parser.pyc and b/pycrs/parser.pyc differ diff --git a/pycrs/projections.py b/pycrs/projections.py index 529bde2..707633d 100644 --- a/pycrs/projections.py +++ b/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" + + + + + + + + + + + + + + + diff --git a/pycrs/projections.pyc b/pycrs/projections.pyc index c6e9fe1..a9a8add 100644 Binary files a/pycrs/projections.pyc and b/pycrs/projections.pyc differ diff --git a/pycrs/units.py b/pycrs/units.py index 6b2a174..6c2cb18 100644 --- a/pycrs/units.py +++ b/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" diff --git a/pycrs/units.pyc b/pycrs/units.pyc index 1220c6e..189d3bf 100644 Binary files a/pycrs/units.pyc and b/pycrs/units.pyc differ diff --git a/tester.py b/tester.py index 0724c3a..28e565d 100644 --- a/tester.py +++ b/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()