diff --git a/pycrs/datums.py b/pycrs/datums.py index 20a1606..0d3cbf4 100644 --- a/pycrs/datums.py +++ b/pycrs/datums.py @@ -1,4 +1,26 @@ + +def find(datumname, crstype, strict=False): + if not strict: + datumname = datumname.lower() + for itemname,item in globals().items(): + if itemname.startswith("_"): + continue + try: + if hasattr(item, crstype): + itemname = getattr(item, crstype) + if not strict: + itemname = itemname.lower() + if datumname == itemname: + return item + except: + pass + else: + return None + + + + class WGS84: proj4 = "WGS84" ogc_wkt = "WGS_1984" @@ -31,6 +53,14 @@ class NAD27: ellipsdef = "" # ellipsoids... to_wgs84 = None +class SphereArcInfo: + proj4 = "" # no name + ogc_wkt = "D_Sphere_ARC_INFO" # confirmed but odd that uses D_ + esri_wkt = "D_Sphere_ARC_INFO" + + semimaj_ax = 6370997.0 + inv_flat = 0.0 + class Unknown: proj4 = "unknown" # no datum name, just ellips + towgs84 params... ogc_wkt = "Unknown" diff --git a/pycrs/datums.pyc b/pycrs/datums.pyc index cdade30..e74ee25 100644 Binary files a/pycrs/datums.pyc and b/pycrs/datums.pyc differ diff --git a/pycrs/ellipsoids.py b/pycrs/ellipsoids.py index 37ad624..da1d133 100644 --- a/pycrs/ellipsoids.py +++ b/pycrs/ellipsoids.py @@ -1,5 +1,27 @@ +def find(ellipsname, crstype, strict=False): + if not strict: + ellipsname = ellipsname.lower() + for itemname,item in globals().items(): + if itemname.startswith("_"): + continue + try: + if hasattr(item, crstype): + itemname = getattr(item, crstype) + if not strict: + itemname = itemname.lower() + if ellipsname == itemname: + return item + except: + pass + else: + return None + + + + + class WGS84: proj4 = "WGS84" ogc_wkt = "WGS_1984" @@ -49,5 +71,20 @@ class Airy1830: semimaj_ax = 6377563.396 inv_flat = 299.3249646 +class SphereArcInfo: + proj4 = "" # no name + ogc_wkt = "Sphere_ARC_INFO" + esri_wkt = "Sphere_ARC_INFO" + + semimaj_ax = 6370997.0 + inv_flat = 0.0 + + + + + + + + diff --git a/pycrs/ellipsoids.pyc b/pycrs/ellipsoids.pyc index 073532d..acd0826 100644 Binary files a/pycrs/ellipsoids.pyc and b/pycrs/ellipsoids.pyc differ diff --git a/pycrs/parameters.py b/pycrs/parameters.py index 25092f4..11f98bc 100644 --- a/pycrs/parameters.py +++ b/pycrs/parameters.py @@ -22,22 +22,44 @@ # +unit and +to_metre are what makes up 'UNIT["Meter",1.0]' +from . import datums from . import directions +################ -# ONLY PURE VALUES INSIDE ELLIPSOID... +def find(paramname, crstype, strict=False): + if not strict: + paramname = paramname.lower() + for itemname,item in globals().items(): + if itemname.startswith("_"): + continue + try: + if hasattr(item, crstype): + itemname = getattr(item, crstype) + if not strict: + itemname = itemname.lower() + if paramname == itemname: + return item + except: + pass + else: + return None + -##+a Semimajor radius of the ellipsoid axis -class SemiMajorRadius: - proj4 = "+a" - def __init__(self, value): - pass +# NOT CURRENTLY USED, BUT SHOULD BE +# FOR NOW ONLY SET AS PURE VALUES INSIDE ELLIPSOID... -##+b Semiminor radius of the ellipsoid axis -class SemiMinorRadius: - proj4 = "+b" - def __init__(self, value): - pass +####+a Semimajor radius of the ellipsoid axis +##class SemiMajorRadius: +## proj4 = "+a" +## def __init__(self, value): +## pass +## +####+b Semiminor radius of the ellipsoid axis +##class SemiMinorRadius: +## proj4 = "+b" +## def __init__(self, value): +## pass @@ -49,6 +71,7 @@ class Azimuth: esri_wkt = "azimuth" ogc_wkt = "azimuth" geotiff = "AzimuthAngle" + def __init__(self, value): self.value = value @@ -59,10 +82,14 @@ class Azimuth: return 'PARAMETER["Azimuth",%s]' % self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["Azimuth",%s]' % self.value ##+datum Datum name (see `proj -ld`) class Datum: + proj4 = "+datum" + ogc_wkt = "DATUM" + esri_wkt = "DATUM" + def __init__(self, name, ellipsoid, datumshift=None): """ Arguments: @@ -77,6 +104,8 @@ class Datum: def to_proj4(self): if self.datumshift: return "%s %s" % (self.ellips.to_proj4(), self.datumshift.to_proj4()) + elif isinstance(self.name, datums.Unknown): + return "%s" % self.ellips.to_proj4() else: return "+datum=%s %s" % (self.name.proj4, self.ellips.to_proj4()) @@ -87,7 +116,10 @@ class Datum: return 'DATUM["%s", %s]' % (self.name.ogc_wkt, self.ellips.to_ogc_wkt()) def to_esri_wkt(self): - return self.to_ogc_wkt() + if self.datumshift: + return 'DATUM["%s", %s, %s]' % (self.name.esri_wkt, self.ellips.to_esri_wkt(), self.datumshift.to_esri_wkt()) + else: + return 'DATUM["%s", %s]' % (self.name.esri_wkt, self.ellips.to_esri_wkt()) def to_geotiff(self): pass @@ -95,6 +127,10 @@ class Datum: ##+ellps Ellipsoid name (see `proj -le`) class Ellipsoid: + proj4 = "+ellps" + ogc_wkt = "SPHEROID" + esri_wkt = "SPHEROID" + def __init__(self, name, semimaj_ax=None, inv_flat=None): """ Arguments: @@ -119,7 +155,7 @@ class Ellipsoid: return 'SPHEROID["%s", %s, %s]' % (self.name.ogc_wkt, self.semimaj_ax, self.inv_flat) def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'SPHEROID["%s", %s, %s]' % (self.name.esri_wkt, self.semimaj_ax, self.inv_flat) def to_geotiff(self): pass @@ -127,6 +163,9 @@ class Ellipsoid: #GEOGCS class GeogCS: + ogc_wkt = "GEOGCS" + esri_wkt = "GEOGCS" + def __init__(self, name, datum, prime_mer, angunit, twin_ax=None): """ Arguments: @@ -143,16 +182,20 @@ class GeogCS: self.twin_ax = twin_ax def to_proj4(self): - return "%s %s %s" % (self.datum.to_proj4(), self.prime_mer.to_proj4(), self.angunit.to_proj4() ) #+axis= AND #, self.twin_ax[0].proj4, self.twin_ax[1].proj4 ) + # dont parse axis to proj4, because in proj4, axis only applies to the cs, ie the projcs (not the geogcs, where wkt can specify with axis) + return "%s %s %s" % (self.datum.to_proj4(), self.prime_mer.to_proj4(), self.angunit.to_proj4() ) def to_ogc_wkt(self): return 'GEOGCS["%s", %s, %s, %s, AXIS["Lon", %s], AXIS["Lat", %s]]' % (self.name, self.datum.to_ogc_wkt(), self.prime_mer.to_ogc_wkt(), self.angunit.to_ogc_wkt(), self.twin_ax[0].ogc_wkt, self.twin_ax[1].ogc_wkt ) def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'GEOGCS["%s", %s, %s, %s, AXIS["Lon", %s], AXIS["Lat", %s]]' % (self.name, self.datum.to_esri_wkt(), self.prime_mer.to_esri_wkt(), self.angunit.to_esri_wkt(), self.twin_ax[0].esri_wkt, self.twin_ax[1].esri_wkt ) #PROJCS class ProjCS: + ogc_wkt = "PROJCS" + esri_wkt = "PROJCS" + def __init__(self, name, geogcs, proj, params, unit, twin_ax=None): """ Arguments: @@ -173,7 +216,6 @@ class ProjCS: 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 @@ -185,11 +227,18 @@ class ProjCS: return string def to_esri_wkt(self): - return self.to_ogc_wkt() + string = 'PROJCS["%s", %s, %s, ' % (self.name, self.geogcs.to_esri_wkt(), self.proj.to_esri_wkt() ) + string += ", ".join(param.to_esri_wkt() for param in self.params) + string += ', %s' % self.unit.to_esri_wkt() + string += ', AXIS["X", %s], AXIS["Y", %s]]' % (self.twin_ax[0].esri_wkt, self.twin_ax[1].esri_wkt ) + return string ##+k Scaling factor (old name) ##+k_0 Scaling factor (new name) class ScalingFactor: + proj4 = "+k" + ogc_wkt = "scale_factor" + def __init__(self, value): self.value = value @@ -200,6 +249,7 @@ class ScalingFactor: return 'PARAMETER["scale_factor", %s]' %self.value def to_esri_wkt(self): + # REALLY?? raise Exception("Paramater not supported by ESRI WKT") def to_geotiff(self): @@ -208,6 +258,10 @@ class ScalingFactor: ##+lat_0 Latitude of origin class LatitudeOrigin: + proj4 = "+lat_0" + ogc_wkt = "latitude_of_origin" + esri_wkt = "Latitude_Of_Origin" + def __init__(self, value): self.value = value @@ -219,7 +273,7 @@ class LatitudeOrigin: return 'PARAMETER["latitude_of_origin", %s]' %self.value def to_esri_wkt(self): - raise Exception("Paramater not supported by ESRI WKT") + return 'PARAMETER["Latitude_Of_Origin", %s]' %self.value def to_geotiff(self): pass @@ -227,6 +281,10 @@ class LatitudeOrigin: ##+lat_1 Latitude of first standard parallel class LatitudeFirstStndParallel: + proj4 = "+lat_1" + ogc_wkt = "standard_parallel_1" + esri_wkt = "Standard_Parallel_1" + def __init__(self, value): self.value = value @@ -237,7 +295,7 @@ class LatitudeFirstStndParallel: return 'PARAMETER["standard_parallel_1", %s]' %self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["Standard_Parallel_1", %s]' %self.value def to_geotiff(self): pass @@ -245,6 +303,10 @@ class LatitudeFirstStndParallel: ##+lat_2 Latitude of second standard parallel class LatitudeSecondStndParallel: + proj4 = "+lat_2" + ogc_wkt = "standard_parallel_2" + esri_wkt = "Standard_Parallel_2" + def __init__(self, value): self.value = value @@ -255,7 +317,7 @@ class LatitudeSecondStndParallel: return 'PARAMETER["standard_parallel_2", %s]' %self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["Standard_Parallel_2", %s]' %self.value def to_geotiff(self): pass @@ -263,6 +325,10 @@ class LatitudeSecondStndParallel: ##+lat_ts Latitude of true scale class LatitudeTrueScale: + proj4 = "lat_ts" + ogc_wkt = "Standard_Parallel_1" + esri_wkt = "Standard_Parallel_1" + def __init__(self, value): self.value = value @@ -273,7 +339,7 @@ class LatitudeTrueScale: return 'PARAMETER["Standard_Parallel_1", %s]' %self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["Standard_Parallel_1", %s]' %self.value def to_geotiff(self): pass @@ -281,6 +347,10 @@ class LatitudeTrueScale: ##+lon_0 Central meridian class CentralMeridian: + proj4 = "+lon_0" + ogc_wkt = "Central_Meridian" + esri_wkt = "Central_Meridian" + def __init__(self, value): self.value = value @@ -291,7 +361,7 @@ class CentralMeridian: return 'PARAMETER["Central_Meridian", %s]' %self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["Central_Meridian", %s]' %self.value def to_geotiff(self): pass @@ -300,6 +370,9 @@ class CentralMeridian: ##+lonc ? Longitude used with Oblique Mercator and possibly a few others class LongitudeCenter: proj4 = "+lonc" + ogc_wkt = "Longitude_Of_Center" + esri_wkt = "Longitude_Of_Center" + def __init__(self, value): self.value = value @@ -310,7 +383,7 @@ class LongitudeCenter: return 'PARAMETER["Longitude_Of_Center", %s]' %self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["Longitude_Of_Center", %s]' %self.value ##+lon_wrap Center longitude to use for wrapping (see below) @@ -318,6 +391,10 @@ class LongitudeCenter: ##+pm Alternate prime meridian (typically a city name, see below) class PrimeMeridian: + proj4 = "+pm" + ogc_wkt = "PRIMEM" + esri_wkt = "PRIMEM" + def __init__(self, value): """ Arguments: @@ -333,10 +410,14 @@ class PrimeMeridian: return 'PRIMEM["Greenwich", %s]' %self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PRIMEM["Greenwich", %s]' %self.value ##+proj Projection name (see `proj -l`) class Projection: + proj4 = "+proj" + ogc_wkt = "PROJECTION" + esri_wkt = "PROJECTION" + def __init__(self, value): self.value = value @@ -347,7 +428,7 @@ class Projection: return 'PROJECTION["%s"]' %self.value.ogc_wkt def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PROJECTION["%s"]' %self.value.esri_wkt ##+zone UTM zone @@ -355,6 +436,9 @@ class Projection: ##+towgs84 3 or 7 term datum transform parameters (see below) class DatumShift: + proj4 = "+towgs84" + ogc_wkt = "TOWGS84" + def __init__(self, value): self.value = value @@ -369,6 +453,8 @@ class DatumShift: ##+to_meter Multiplier to convert map units to 1.0m class MeterMultiplier: + proj4 = "+to_meter" + def __init__(self, value): self.value = value @@ -380,10 +466,12 @@ class MeterMultiplier: return bytes(self.value) def to_esri_wkt(self): - return self.to_ogc_wkt() + return bytes(self.value) ##+units meters, US survey feet, etc. class UnitType: + proj4 = "+units" + def __init__(self, value): """ Arguments: @@ -400,10 +488,13 @@ class UnitType: return bytes(self.value.ogc_wkt) def to_esri_wkt(self): - return self.to_ogc_wkt() + return bytes(self.value.esri_wkt) # special... class Unit: + ogc_wkt = "UNIT" + esri_wkt = "UNIT" + def __init__(self, unittype, metermultiplier): self.unittype = unittype self.metermultiplier = metermultiplier @@ -415,10 +506,13 @@ class Unit: return 'UNIT["%s", %s]' %(self.unittype.to_ogc_wkt(), self.metermultiplier.to_ogc_wkt()) def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'UNIT["%s", %s]' %(self.unittype.to_esri_wkt(), self.metermultiplier.to_esri_wkt()) # angular unit class AngularUnit: + ogc_wkt = "UNIT" + esri_wkt = "UNIT" + def __init__(self, unittype, metermultiplier): self.unittype = unittype self.metermultiplier = metermultiplier @@ -431,7 +525,7 @@ class AngularUnit: return 'UNIT["%s", %s]' %(self.unittype.to_ogc_wkt(), self.metermultiplier.to_ogc_wkt()) def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'UNIT["%s", %s]' %(self.unittype.to_esri_wkt(), self.metermultiplier.to_esri_wkt()) ##+x_0 False easting class FalseEasting: @@ -450,7 +544,7 @@ class FalseEasting: return 'PARAMETER["false_easting", %s]' % self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["False_Easting", %s]' % self.value ##+y_0 False northing class FalseNorthing: @@ -469,10 +563,14 @@ class FalseNorthing: return 'PARAMETER["false_northing", %s]' % self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["False_Northing", %s]' % self.value ##+h Satellite height class SatelliteHeight: + proj4 = "+h" + ogc_wkt = "satellite_height" + esri_wkt = "satellite_height" + def __init__(self, value): self.value = value @@ -483,10 +581,12 @@ class SatelliteHeight: return 'PARAMETER["satellite_height", %s]' % self.value def to_esri_wkt(self): - return self.to_ogc_wkt() + return 'PARAMETER["satellite_height", %s]' % self.value ##+tilt Tilt angle class TiltAngle: + proj4 = "+tilt" + def __init__(self, value): self.value = value @@ -494,10 +594,10 @@ class TiltAngle: return "+tilt=%s" %self.value def to_ogc_wkt(self): - raise Exception("Paramater not supported by OGC WKT") + raise Exception("Parameter not supported by OGC WKT") def to_esri_wkt(self): - raise Exception("Paramater not supported by ESRI WKT") + raise Exception("Parameter not supported by ESRI WKT") # then the final CRS object which is instantiated with all of these? @@ -508,12 +608,15 @@ class CRS: self.toplevel = toplevel def to_proj4(self): - return "%s +no_defs" % self.toplevel.to_proj4() + if isinstance(self.toplevel, ProjCS): + return "%s +no_defs" % self.toplevel.to_proj4() + elif isinstance(self.toplevel, GeogCS): + return "+proj=longlat %s +no_defs" % self.toplevel.to_proj4() def to_ogc_wkt(self): return "%s" % self.toplevel.to_ogc_wkt() def to_esri_wkt(self): - return self.to_ogc_wkt() + return "%s" % self.toplevel.to_esri_wkt() diff --git a/pycrs/parameters.pyc b/pycrs/parameters.pyc index df59ae7..9da5b5e 100644 Binary files a/pycrs/parameters.pyc and b/pycrs/parameters.pyc differ diff --git a/pycrs/parser.pyc b/pycrs/parser.pyc index 010afc8..4eb4c9b 100644 Binary files a/pycrs/parser.pyc and b/pycrs/parser.pyc differ diff --git a/pycrs/projections.py b/pycrs/projections.py index 707633d..a24c1cb 100644 --- a/pycrs/projections.py +++ b/pycrs/projections.py @@ -1,20 +1,46 @@ +def find(projname, crstype, strict=False): + if not strict: + projname = projname.lower() + for itemname,item in globals().items(): + if itemname.startswith("_"): + continue + try: + if hasattr(item, crstype): + itemname = getattr(item, crstype) + if not strict: + itemname = itemname.lower() + if projname == itemname: + return item + except: + pass + else: + return None + + + + + class Robinson: proj4 = "robin" ogc_wkt = "Robinson" + esri_wkt = "Robinson" class UTM: proj4 = "utm" ogc_wkt = "Transverse_Mercator" + esri_wkt = "Transverse_Mercator" class ObliqueMercator: proj4 = "omerc" ogc_wkt = "Hotine_Oblique_Mercator" + esri_wkt = "Hotine_Oblique_Mercator_Azimuth_Natural_Origin" class AlbersEqualArea: proj4 = "aea" - ogc_wkt = "Albers" + ogc_wkt = "Albers_Conic_Equal_Area" + esri_wkt = "Albers" class CylindricalEqualArea: proj4 = "cea" @@ -24,90 +50,112 @@ class CylindricalEqualArea: class EquiDistantConic: proj4 = "eqdc" ogc_wkt = "Equidistant_Conic" + esri_wkt = "Equidistant_Conic" class TransverseMercator: proj4 = "tmerc" ogc_wkt = "Transverse_Mercator" + esri_wkt = "Transverse_Mercator" class GallStereographic: proj4 = "gall" ogc_wkt = "Gall_Stereographic" + esri_wkt = "Gall_Stereographic" class Gnomonic: proj4 = "gnom" ogc_wkt = "Gnomonic" + esri_wkt = "Gnomonic" class LambertAzimuthalEqualArea: proj4 = "laea" ogc_wkt = "Lambert_Azimuthal_Equal_Area" + esri_wkt = "Lambert_Azimuthal_Equal_Area" class MillerCylindrical: proj4 = "mill" ogc_wkt = "Miller_Cylindrical" + esri_wkt = "Miller_Cylindrical" class Mollweide: proj4 = "moll" ogc_wkt = "Mollweide" + esri_wkt = "Mollweide" class ObliqueStereographic: proj4 = "sterea" ogc_wkt = "Oblique_Stereographic" + esri_wkt = "Stereographic_North_Pole" class Orthographic: proj4 = "ortho" ogc_wkt = "Orthographic" + esri_wkt = "Orthographic" class PolarStereographic: proj4 = "stere" - ogc_wkt = "Polar_Stereographic" + ogc_wkt = "Polar_Stereographic" # could also be just stereographic + esri_wkt = "Stereographic" # but also spelled with additional _South/North_Pole, for the same projection and diff params (maybe just for humans)?... class Sinusoidal: proj4 = "sinu" ogc_wkt = "Sinusoidal" + esri_wkt = "Sinusoidal" class VanDerGrinten: proj4 = "vandg" ogc_wkt = "VanDerGrinten" + esri_wkt = "Van_der_Grinten_I" class Equirectangular: proj4 = "eqc" ogc_wkt = "Equirectangular" + esri_wkt = "Equidistant_Cylindrical" class LambertConformalConic: proj4 = "lcc" - ogc_wkt = "Lambert_Conformal_Conic" + ogc_wkt = "Lambert_Conformal_Conic" # possible has some variants + esri_wkt = "Lambert_Conformal_Conic" class Krovak: proj4 = "krovak" ogc_wkt = "Krovak" + esri_wkt = "Krovak" class NearSidedPerspective: proj4 = "nsper" ogc_wkt = "Near_sided_perspective" + esri_wkt = "Near_sided_perspective" # not confirmed class TiltedPerspective: proj4 = "tsper" ogc_wkt = "Tilted_perspective" + esri_wkt = "Tilted_perspective" # not confirmed class InteruptedGoodeHomolosine: proj4 = "igh" ogc_wkt = "Interrupted_Goodes_Homolosine" + esri_wkt = "Interrupted_Goodes_Homolosine" class Larrivee: proj4 = "larr" ogc_wkt = "Larrivee" + esri_wkt = "Larrivee" # not confirmed class LamberEqualAreaConic: proj4 = "leac" ogc_wkt = "Lambert_Equal_Area_Conic" + esri_wkt = "Lambert_Equal_Area_Conic" # not confirmed class Mercator: proj4 = "merc" - ogc_wkt = "Mercator" + ogc_wkt = "Mercator" # has multiple varieties + esri_wkt = "Mercator" class ObliqueCylindricalEqualArea: proj4 = "ocea" ogc_wkt = "Oblique_Cylindrical_Equal_Area" + esri_wkt = "Oblique_Cylindrical_Equal_Area" diff --git a/pycrs/projections.pyc b/pycrs/projections.pyc index a9a8add..8b7c509 100644 Binary files a/pycrs/projections.pyc and b/pycrs/projections.pyc differ diff --git a/pycrs/units.py b/pycrs/units.py index 6c2cb18..9b545f7 100644 --- a/pycrs/units.py +++ b/pycrs/units.py @@ -1,5 +1,30 @@ +def find(unitname, crstype, strict=False): + if not strict: + unitname = unitname.lower() + for itemname,item in globals().items(): + if itemname.startswith("_"): + continue + try: + if hasattr(item, crstype): + itemname = getattr(item, crstype) + # compare + if not strict: + itemname = itemname.lower() + if unitname == itemname: + return item + # special handling of wkt meters which has multiple possibilities + elif isinstance(itemname, Meter) and crstype.endswith("wkt") and not strict and unitname in ("meters","meter","metre"): + return item + except: + pass + else: + return None + + + + class Meter: proj4 = "m" ogc_wkt = "Meters" # or is it metre?? sometimes even Meter? diff --git a/pycrs/units.pyc b/pycrs/units.pyc index 189d3bf..814c2f7 100644 Binary files a/pycrs/units.pyc and b/pycrs/units.pyc differ diff --git a/pycrs/webscrape.pyc b/pycrs/webscrape.pyc index ca30c96..dd753cd 100644 Binary files a/pycrs/webscrape.pyc and b/pycrs/webscrape.pyc differ diff --git a/tester.py b/tester.py index e3e8e43..36f3bf1 100644 --- a/tester.py +++ b/tester.py @@ -12,16 +12,15 @@ def render_world(crs): import pyagg import pyproj import random - c = pyagg.Canvas(600,600) - c.geographic_space() + + # load world borders 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 + # convert coordinates + fromproj = pyproj.Proj("+init=EPSG:4326") + toproj = pyproj.Proj(crs.to_proj4()) for feat in data: if feat.geometry.type == "Polygon": feat.geometry.coordinates = [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1])) @@ -31,19 +30,33 @@ def render_world(crs): [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1])) for ring in poly] for poly in feat.geometry.coordinates] + # get zoom area + # NOTE: PYGEOJ data.bbox doesnt update.........FIX! 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 + data.update_bbox() + bbox = data.bbox + +## # to avoid inf bounds and no render in satellite view +## xmins, ymins, xmaxs, ymaxs = zip(*(feat.geometry.bbox for feat in data)) +## inf = float("inf") +## xmaxs = (xmax for xmax in xmaxs if xmax != inf) +## ymaxs = (ymax for ymax in ymaxs if ymax != inf) +## bbox = (min(xmins), min(ymins), max(xmaxs), max(ymaxs)) + + # set up drawing + c = pyagg.Canvas(600,600) + c.geographic_space() c.zoom_bbox(*bbox) + c.zoom_out(1.3) + + # draw 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 + print("unable to draw?", feat.geometry) c.view() @@ -53,27 +66,41 @@ def render_world(crs): ########################### # OGC WKT -print "testing ogc wkt" -crs = pycrs.parser.from_sr_code(54030) +print("--------") +print("testing ogc wkt") +print("") +#crs = pycrs.parser.from_sr_code(54030) #crs = pycrs.parser.from_sr_code(22) -print crs.to_proj4() +#crs = pycrs.parser.from_esri_code(54031) +crs = pycrs.parser.from_sr_code(6978) wkt = crs.to_ogc_wkt() -print "Original:", wkt -print "Reconstructed:", pycrs.parser.from_ogc_wkt(wkt).to_ogc_wkt() -render_world(crs) +print("Original:", wkt) +print("") + +crs = pycrs.parser.from_ogc_wkt(wkt) + +print("Reconstructed:", crs.to_ogc_wkt()) +print("") + +#render_world(crs) ########################### # PROJ4 -print "testing proj4" -#proj4 = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" +print("--------") +print("testing proj4") +print("") +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=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=-60 +lat_0=40 +h=2000000000000000000000000" + +print("Original:", proj4) +print("") #crs = pycrs.parser.from_esri_code(54030) #crs = pycrs.parser.from_sr_code(7898) @@ -82,17 +109,31 @@ proj4 = "+proj=larr +datum=WGS84 +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 crs = pycrs.parser.from_proj4(proj4) #crs = pycrs.parser.from_ogc_wkt(wkt) -print crs.to_proj4() +print("Reconstructed:", crs.to_proj4()) +print("") -render_world(crs) +#render_world(crs) ########################### -# PRJ FILE +# ESRI WKT/PRJ FILE +print("--------") +print("testing esri prj file") +print("") crs = pycrs.loader.from_file("testfiles/natearth.prj") -print crs +wkt = crs.to_esri_wkt() + +print("Original:", wkt) +print("") + +crs = pycrs.parser.from_esri_wkt(wkt) + +print("Reconstructed:", crs.to_esri_wkt()) +print("") + +#render_world(crs)