Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- lyr = conn.GetLayerByName(tbl) # Where conn is OGR PG connection
- srs = ly.GetSpatialRef()
- print srs
- PROJCS["OSGB 1936 / British National Grid",
- GEOGCS["OSGB 1936",
- DATUM["OSGB_1936",
- SPHEROID["Airy 1830",6377563.396,299.3249646,
- AUTHORITY["EPSG","7001"]],
- AUTHORITY["EPSG","6277"]],
- PRIMEM["Greenwich",0,
- AUTHORITY["EPSG","8901"]],
- UNIT["degree",0.01745329251994328,
- AUTHORITY["EPSG","9122"]],
- AUTHORITY["EPSG","4277"]],
- UNIT["metre",1,
- AUTHORITY["EPSG","9001"]],
- PROJECTION["Transverse_Mercator"],
- PARAMETER["latitude_of_origin",49],
- PARAMETER["central_meridian",-2],
- PARAMETER["scale_factor",0.9996012717],
- PARAMETER["false_easting",400000],
- PARAMETER["false_northing",-100000],
- AUTHORITY["EPSG","27700"],
- AXIS["Easting",EAST],
- AXIS["Northing",NORTH]]
- srs.GetEPSG()
- print srs
- 27700
- In [1]: import osgeo.osr as osr
- In [2]: srs = osr.SpatialReference()
- In [3]: srs.SetFromUserInput("EPSG:27700")
- Out[3]: 0
- In [4]: print srs
- PROJCS["OSGB 1936 / British National Grid",
- GEOGCS["OSGB 1936",
- DATUM["OSGB_1936",
- SPHEROID["Airy 1830",6377563.396,299.3249646,
- AUTHORITY["EPSG","7001"]],
- TOWGS84[375,-111,431,0,0,0,0],
- AUTHORITY["EPSG","6277"]],
- PRIMEM["Greenwich",0,
- AUTHORITY["EPSG","8901"]],
- UNIT["degree",0.0174532925199433,
- AUTHORITY["EPSG","9122"]],
- AUTHORITY["EPSG","4277"]],
- PROJECTION["Transverse_Mercator"],
- PARAMETER["latitude_of_origin",49],
- PARAMETER["central_meridian",-2],
- PARAMETER["scale_factor",0.9996012717],
- PARAMETER["false_easting",400000],
- PARAMETER["false_northing",-100000],
- UNIT["metre",1,
- AUTHORITY["EPSG","9001"]],
- AXIS["Easting",EAST],
- AXIS["Northing",NORTH],
- AUTHORITY["EPSG","27700"]]
- In [5]: srs.GetAttrValue("AUTHORITY", 0)
- Out[5]: 'EPSG'
- In [6]: srs.GetAttrValue("AUTHORITY", 1)
- Out[6]: '27700'
- In [12]: srs.GetAttrValue("PRIMEM|AUTHORITY", 1)
- Out[12]: '8901'
- In [13]: srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
- Out[13]: '4277'
- def wkt2epsg(wkt, epsg='/usr/local/share/proj/epsg', forceProj4=False):
- ''' Transform a WKT string to an EPSG code
- Arguments
- ---------
- wkt: WKT definition
- epsg: the proj.4 epsg file (defaults to '/usr/local/share/proj/epsg')
- forceProj4: whether to perform brute force proj4 epsg file check (last resort)
- Returns: EPSG code
- '''
- code = None
- p_in = osr.SpatialReference()
- s = p_in.ImportFromWkt(wkt)
- if s == 5: # invalid WKT
- return None
- if p_in.IsLocal() == 1: # this is a local definition
- return p_in.ExportToWkt()
- if p_in.IsGeographic() == 1: # this is a geographic srs
- cstype = 'GEOGCS'
- else: # this is a projected srs
- cstype = 'PROJCS'
- an = p_in.GetAuthorityName(cstype)
- ac = p_in.GetAuthorityCode(cstype)
- if an is not None and ac is not None: # return the EPSG code
- return '%s:%s' %
- (p_in.GetAuthorityName(cstype), p_in.GetAuthorityCode(cstype))
- else: # try brute force approach by grokking proj epsg definition file
- p_out = p_in.ExportToProj4()
- if p_out:
- if forceProj4 is True:
- return p_out
- f = open(epsg)
- for line in f:
- if line.find(p_out) != -1:
- m = re.search('<(\d+)>', line)
- if m:
- code = m.group(1)
- break
- if code: # match
- return 'EPSG:%s' % code
- else: # no match
- return None
- else:
- return None
- In [1]: import osgeo.osr as osr
- In [2]: srs = osr.SpatialReference()
- In [3]: srs.SetFromUserInput("EPSG:27700")
- Out[3]: 0
- In [4]: srs.GetAuthorityCode(None)
- Out[4]: '27700'
- In [5]: srs.GetAuthorityCode(None)
- Out[5]: 'EPSG'
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement