Basic fromproj4 works, except need better handle when unprojected (longlat and latlong)

Also just need to expand to more projection name, datum name, and
ellipsoid name definitions, and also additional special parameters.

After that, test, and move onto from_wkt()...
This commit is contained in:
Karim Bahgat
2015-06-29 03:29:02 +02:00
parent 5a46702471
commit fcc463fa76
24 changed files with 504 additions and 132 deletions
-12
View File
@@ -1,12 +0,0 @@
# CRS Names
class proj_NAD87:
proj4 = ""
ogc_wkt = ""
class datum_WGS84:
proj4 = ""
ogc_wkt = ""
-88
View File
@@ -1,88 +0,0 @@
# parse from text strings
# 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
def from_epsg_code(string):
# must go online (or look up local table) to get crs details
pass
def from_esri_code(string):
# must go online (or look up local table) to get crs details
pass
def from_sr_code(string):
# must go online (or look up local table) to get crs details
pass
def from_esri_wkt(string):
# 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_unknown_wkt(string):
# detect if ogc wkt or esri wkt
# TIPS: esri wkt datums all use "D_" before the datum name
# then load with appropriate function
pass
def from_proj4(string):
# parse arguments into components
# use args to create crs
pass
def from_ogc_urn(string):
# hmmm, seems like ogc urn could be anything incl online link, epsg, etc...
# if necessary, must go online (or lookup local table) to get details
# maybe test which of these and run their function?
# examples urn:ogc:def:crs:OGC:1.3:CRS1
# or with EPSG instead of OGC
# If OGC, 1.3 is pdf version, and after that is a name from list below
# as found in pdf: "Definition identifier URNs in OGC namespace"
# OGC crs definitions
# URN | CRS name | Definition reference
# urn:ogc:def:crs:OGC:1.3:CRS1 Map CS B.2 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS84 WGS 84 longitude-latitude B.3 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS83 NAD83 longitude-latitude B.4 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS27 NAD27 longitude-latitude B.5 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS88 NAVD 88 B.6 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42001:99:8888 Auto universal transverse mercator B.7 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42002:99:8888 Auto transverse mercator B.8 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42003:99:8888 Auto orthographic B.9 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42004:99:8888 Auto equirectangular B.10 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42005:99 Auto Mollweide B.11 in OGC 06-042
pass
def from_unknown_text(text):
# detect type and load with appropriate function
if string.startswith("urn:"):
from_ogc_urn(string)
elif string.startswith("+proj="):
from_proj4(string)
elif string.startswith("PROJCS["):
from_unknown_wkt(string)
elif string.startswith("EPSG:"):
from_epsg_code(string)
elif string.startswith("ESRI:"):
from_esri_code(string)
elif string.startswith("SR-ORG:"):
from_sr_code(string)
else: raise Exception("Could not detect which type of crs")
def from_geotiff_params(**params):
+4
View File
@@ -0,0 +1,4 @@
from . import loader
from . import parser
from . import webscrape
Binary file not shown.
View File
+5
View File
@@ -0,0 +1,5 @@
class WGS84:
proj4 = "WGS84"
ogc_wkt = "WGS_1984"
BIN
View File
Binary file not shown.
+31
View File
@@ -0,0 +1,31 @@
class North:
proj4 = "n"
ogc_wkt = "NORTH"
esri_wkt = "NORTH"
class East:
proj4 = "e"
ogc_wkt = "EAST"
esri_wkt = "EAST"
class South:
proj4 = "s"
ogc_wkt = "SOUTH"
esri_wkt = "SOUTH"
class West:
proj4 = "w"
ogc_wkt = "WEST"
esri_wkt = "WEST"
class Up:
proj4 = "u"
ogc_wkt = "UP"
esri_wkt = "UP"
class Down:
proj4 = "d"
ogc_wkt = "DOWN"
esri_wkt = "DOWN"
Binary file not shown.
+9
View File
@@ -0,0 +1,9 @@
class WGS84:
proj4 = "WGS84"
ogc_wkt = "WGS_1984"
semimaj_ax = 6378137
inv_flat = 298.257223563
Binary file not shown.
View File
BIN
View File
Binary file not shown.
+164 -32
View File
@@ -22,6 +22,8 @@
# +unit and +to_metre are what makes up 'UNIT["Meter",1.0]'
from . import directions
##+a Semimajor radius of the ellipsoid axis
@@ -39,15 +41,6 @@ class Azimuth:
def __init__(self, value):
pass
##+axis Axis orientation (new in 4.8.0)
class AxisOrientation:
proj4 = "+axis"
esri_wkt = None
ogc_wkt = "AXIS"
geotiff = None
def __init__(self, value):
pass
##+b Semiminor radius of the ellipsoid axis
class SemiMinorRadius:
proj4 = "+b"
@@ -55,15 +48,22 @@ class SemiMinorRadius:
pass
##+datum Datum name (see `proj -ld`)
class DatumName:
def __init__(self, value):
self.value = value
class Datum:
def __init__(self, name, ellipsoid):
"""
Arguments:
- **name**: Specific datum name instance.
- **ellipsoid**: Ellipsoid parameter instance.
"""
self.name = name
self.ellips = ellipsoid
def to_proj4(self):
return "+datum=%s" %self.value
return "+datum=%s %s" % (self.name.proj4, self.ellips.to_proj4())
def to_ogc_wkt(self):
return 'DATUM[%s]' %self.value
return 'DATUM["%s", %s]' % (self.name.ogc_wkt, self.ellips.to_ogc_wkt())
def to_esri_wkt(self):
return self.to_ogc_wkt()
@@ -73,22 +73,96 @@ class DatumName:
#return "GeogGeodeticDatum"
##+ellps Ellipsoid name (see `proj -le`)
class EllipsoidName:
def __init__(self, value):
self.value = value
class Ellipsoid:
def __init__(self, name, semimaj_ax=None, inv_flat=None):
"""
Arguments:
- **name**: Specific ellipsoid name instance.
"""
self.name = name
# get default values if not specified
if semimaj_ax == None:
semimaj_ax = self.name.semimaj_ax
if inv_flat == None:
inv_flat = self.name.inv_flat
self.semimaj_ax = semimaj_ax
self.inv_flat = inv_flat
def to_proj4(self):
return "+ellps=%s" %self.value
return "+ellps=%s +a=%s +f=%s" % (self.name.proj4, self.semimaj_ax, self.inv_flat)
def to_ogc_wkt(self):
return 'SPHEROID[%s]' %self.value
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()
def to_geotiff(self):
pass
#return "GeogEllipsoid"
#return "GeogEllipsoid"
#GEOGCS
class GeogCS:
def __init__(self, name, datum, prime_mer, angunit, twin_ax=None):
"""
Arguments:
- **name**: Arbitrary name.
"""
self.name = name
self.datum = datum
self.prime_mer = prime_mer
self.angunit = angunit
if twin_ax == None:
# default axes
twin_ax = directions.East(), directions.North()
self.twin_ax = twin_ax
def to_proj4(self):
# axis excluded because not sure if should be set from geogcs or projcs
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 )
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()
#PROJCS
class ProjCS:
def __init__(self, name, geogcs, proj, params, unit, twin_ax=None):
"""
Arguments:
- **name**: Arbitrary name.
"""
self.name = name
self.geogcs = geogcs
self.proj = proj
self.params = params
if twin_ax == None:
# default axes
twin_ax = directions.East(), directions.North()
self.twin_ax = twin_ax
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)
# axis excluded because not sure if should be set from geogcs or projcs
#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
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 += ', AXIS["X", %s], AXIS["Y", %s]]' % (self.twin_ax[0].ogc_wkt, self.twin_ax[1].ogc_wkt )
return string
def to_esri_wkt(self):
return self.to_ogc_wkt()
##+k Scaling factor (old name)
@@ -190,6 +264,11 @@ class LongitudeSpecial:
##+pm Alternate prime meridian (typically a city name, see below)
class PrimeMeridian:
def __init__(self, value):
"""
Arguments:
- **value**: Longitude value relative to Greenwich.
"""
self.value = value
def to_proj4(self):
@@ -202,15 +281,15 @@ class PrimeMeridian:
return self.to_ogc_wkt()
##+proj Projection name (see `proj -l`)
class ProjectionName:
class Projection:
def __init__(self, value):
self.value = value
def to_proj4(self):
return "+proj=%s" %self.value
return "+proj=%s" %self.value.proj4
def to_ogc_wkt(self):
return 'PROJECTION["%s"]' %self.value
return 'PROJECTION["%s"]' %self.value.ogc_wkt
def to_esri_wkt(self):
return self.to_ogc_wkt()
@@ -250,14 +329,19 @@ class MeterMultiplier:
##+units meters, US survey feet, etc.
class UnitType:
def __init__(self, value):
"""
Arguments:
- **value**: A specific unit type instance, eg Meter().
"""
self.value = value
def to_proj4(self):
return "+units=%s" %self.value
return "+units=%s" %self.value.proj4
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)
return str(self.value.ogc_wkt)
def to_esri_wkt(self):
return self.to_ogc_wkt()
@@ -269,7 +353,23 @@ class Unit:
self.metermultiplier = metermultiplier
def to_proj4(self):
return "+%s +%s" %(self.unittype, self.metermultiplier)
return "%s %s" %(self.unittype.to_proj4(), self.metermultiplier.to_proj4())
def to_ogc_wkt(self):
return 'UNIT["%s", %s]' %(self.unittype.to_ogc_wkt(), self.metermultiplier.to_ogc_wkt())
def to_esri_wkt(self):
return self.to_ogc_wkt()
# angular unit
class AngularUnit:
def __init__(self, unittype, metermultiplier):
self.unittype = unittype
self.metermultiplier = metermultiplier
def to_proj4(self):
# cannot be specified in proj4, so just return nothing
return ""
def to_ogc_wkt(self):
return 'UNIT["%s", %s]' %(self.unittype.to_ogc_wkt(), self.metermultiplier.to_ogc_wkt())
@@ -283,23 +383,55 @@ class FalseEasting:
esri_wkt = "False_Easting"
ogc_wkt = "false_easting"
geotiff = "FalseEasting"
def __init__(self, value):
pass
self.value = value
def to_proj4(self):
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):
return self.to_ogc_wkt()
##+y_0 False northing
class FalseEasting:
class FalseNorthing:
proj4 = "+y_0"
esri_wkt = "False_Northing"
ogc_wkt = "false_northing"
geotiff = "FalseNorthing"
def __init__(self, value):
pass
self.value = value
def to_proj4(self):
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()
# then the final CRS object which is instantiated with all of these?
# remember to use +no_defs when outputting to proj4
# ...
class CRS:
def __init__(self, toplevel):
self.toplevel = toplevel
def to_proj4(self):
return "%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()
Binary file not shown.
+260
View File
@@ -0,0 +1,260 @@
# parse from text strings
# 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
from . import datums
from . import ellipsoids
from . import parameters
from . import units
from . import projections
def from_epsg_code(string):
# must go online (or look up local table) to get crs details
pass
def from_esri_code(string):
# must go online (or look up local table) to get crs details
pass
def from_sr_code(string):
# must go online (or look up local table) to get crs details
pass
def from_esri_wkt(string):
# 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_unknown_wkt(string):
# detect if ogc wkt or esri wkt
# TIPS: esri wkt datums all use "D_" before the datum name
# then load with appropriate function
pass
def from_proj4(string):
# parse arguments into components
# use args to create crs
## proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
##
## PROJCS["World_Robinson",
## GEOGCS["GCS_WGS_1984",
## DATUM["WGS_1984",
## SPHEROID["WGS_1984",6378137,298.257223563]],
## PRIMEM["Greenwich",0],
## UNIT["Degree",0.017453292519943295]],
## PROJECTION["Robinson"],
## PARAMETER["False_Easting",0],
## PARAMETER["False_Northing",0],
## PARAMETER["Central_Meridian",0],
## UNIT["Meter",1],
## AUTHORITY["EPSG","54030"]]
params = []
partdict = dict([part.split("=") for part in string.split()
if not part.startswith("+no_defs")])
# DATUM
# datum param is required
if "+datum" in partdict:
# get predefined datum def
if partdict["+datum"] == "WGS84":
datumdef = datums.WGS84()
# ELLIPS
# ellipse param is required
if "+ellps" in partdict:
# get predefined ellips def
if partdict["+ellps"] == "WGS84":
ellipsdef = ellipsoids.WGS84()
else:
raise Exception("Could not find required +ellps element")
## create datum and ellips param objs
ellips = parameters.Ellipsoid(ellipsdef,
semimaj_ax=partdict.get("+a"),
inv_flat=partdict.get("+f"))
datum = parameters.Datum(datumdef, ellips)
else:
raise Exception("Could not find required +datum element")
# PRIME MERIDIAN
# set default
prime_mer = parameters.PrimeMeridian(0)
# overwrite with user input
if "+pm" in partdict:
# for now only support longitude, later add name support:
## greenwich 0dE
## lisbon 9d07'54.862"W
## paris 2d20'14.025"E
## bogota 74d04'51.3"E
## madrid 3d41'16.48"W
## rome 12d27'8.4"E
## bern 7d26'22.5"E
## jakarta 106d48'27.79"E
## ferro 17d40'W
## brussels 4d22'4.71"E
## stockholm 18d3'29.8"E
## athens 23d42'58.815"E
## oslo 10d43'22.5"E
prime_mer = parameters.PrimeMeridian(partdict["+pm"])
# ANGULAR UNIT
## proj4 cannot set angular unit, so just set to default
metmulti = parameters.MeterMultiplier(0.017453292519943295)
unittype = parameters.UnitType(units.Degree())
angunit = parameters.AngularUnit(unittype, metmulti)
# GEOGCS (note, currently does not load axes)
geogcs = parameters.GeogCS("Unknown", datum, prime_mer, angunit) #, twin_ax)
# PROJECTION
if "+proj" in partdict:
# get predefined proj def
if partdict["+proj"] == "robin":
projdef = projections.Robinson()
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...
# ...
else:
projdef = None
else:
raise Exception("Could not find required +proj element")
if projdef:
# create proj param obj
proj = parameters.Projection(projdef)
# CENTRAL MERIDIAN
if "+lon_0" in partdict:
val = partdict["+lon_0"]
obj = parameters.CentralMeridian(val)
params.append(obj)
# FALSE EASTING
if "+x_0" in partdict:
val = partdict["+x_0"]
obj = parameters.FalseEasting(val)
params.append(obj)
# FALSE NORTHING
if "+y_0" in partdict:
val = partdict["+y_0"]
obj = parameters.FalseNorthing(val)
params.append(obj)
# UNIT
## set default
metmulti = parameters.MeterMultiplier(1.0)
unittype = parameters.UnitType(units.Meter())
## override with user input
if "+to_meter" in partdict:
metmulti = parameters.MeterMultiplier(partdict["+to_meter"])
if "+units" in partdict:
if partdict["+units"] == "m":
unittype = parameters.UnitType(units.Meter())
## create unitobj
unit = parameters.Unit(unittype, metmulti)
# PROJCS
projcs = parameters.ProjCS("Unknown", geogcs, proj, params, unit)
# CRS
crs = parameters.CRS(projcs)
else:
crs = parameters.CRS(geogcs)
# FINISHED
return crs
def from_ogc_urn(string):
# hmmm, seems like ogc urn could be anything incl online link, epsg, etc...
# if necessary, must go online (or lookup local table) to get details
# maybe test which of these and run their function?
# examples urn:ogc:def:crs:OGC:1.3:CRS1
# or with EPSG instead of OGC
# If OGC, 1.3 is pdf version, and after that is a name from list below
# as found in pdf: "Definition identifier URNs in OGC namespace"
# OGC crs definitions
# URN | CRS name | Definition reference
# urn:ogc:def:crs:OGC:1.3:CRS1 Map CS B.2 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS84 WGS 84 longitude-latitude B.3 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS83 NAD83 longitude-latitude B.4 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS27 NAD27 longitude-latitude B.5 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:CRS88 NAVD 88 B.6 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42001:99:8888 Auto universal transverse mercator B.7 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42002:99:8888 Auto transverse mercator B.8 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42003:99:8888 Auto orthographic B.9 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42004:99:8888 Auto equirectangular B.10 in OGC 06-042
# urn:ogc:def:crs:OGC:1.3:AUTO42005:99 Auto Mollweide B.11 in OGC 06-042
pass
def from_unknown_text(text):
# detect type and load with appropriate function
if string.startswith("urn:"):
from_ogc_urn(string)
elif string.startswith("+proj="):
from_proj4(string)
elif string.startswith("PROJCS["):
from_unknown_wkt(string)
elif string.startswith("EPSG:"):
from_epsg_code(string)
elif string.startswith("ESRI:"):
from_esri_code(string)
elif string.startswith("SR-ORG:"):
from_sr_code(string)
else: raise Exception("Could not detect which type of crs")
def from_geotiff_parameters(**params):
pass
BIN
View File
Binary file not shown.
+10
View File
@@ -0,0 +1,10 @@
class Robinson:
proj4 = "robin"
ogc_wkt = "Robinson"
class UTM:
proj4 = "utm"
ogc_wkt = "Transverse_Mercator"
Binary file not shown.
+11
View File
@@ -0,0 +1,11 @@
class Meter:
proj4 = "m"
ogc_wkt = "Meters" # or is it metre?? sometimes even Meter?
esri_wkt = "Meter"
class Degree:
proj4 = "degrees"
ogc_wkt = "degree"
esri_wkt = "Degree"
BIN
View File
Binary file not shown.
View File
Binary file not shown.
+10
View File
@@ -0,0 +1,10 @@
import pycrs
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"
crs = pycrs.parser.from_proj4(proj4)
print crs.to_ogc_wkt()