例如:
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]]
那么我如何获得投影的EPSG值?例如:
srs.GetEPSG()
print srs
27700
我已经尝试过
srs.GetAttrValue('AUTHORITY')
,但这只会返回'EPSG'
。#1 楼
它有点埋头,但是GetAttrValue()有第二个参数,该参数返回该序数的值。因此,我可以这样做: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'
哪些可用于查找投影CS的地理坐标系:
In [13]: srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
Out[13]: '4277'
评论
谢谢,太好了。我将执行它。我没有时间进行进一步的“玩耍”了-缺少GDAL / OGR文档再次阻碍了快速的应用程序开发!
– Tomas
2012年2月14日14:14
我尝试使用带有“ AUTHORITY”和“ 1”参数的GetAttrValue函数,并注意到它并不总是返回EPSG代码,因为EPSG代码并不总是包含在WKT中。对于这种情况,我还是有点模糊。我发现以下解决方案可以很好地满足我的需求:gis.stackexchange.com/questions/7608/…
–挖洞
17年6月28日在20:42
#2 楼
这是一个对我有用的代码片段: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
#3 楼
SpatialReference.GetAuthorityCode()
将None
作为参数,该参数在根元素上找到一个权限节点(即,适当时是投影的/地理的)。同样适用于GetAuthorityName()
: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.GetAuthorityName(None)
Out[5]: 'EPSG'
评论
None表示它将自行在根目录中搜索以找到授权码(在此情况下为epsg)
–尼克
20年4月6日在12:37
评论
我已经尝试过srs.GetAttrValue('AUTHORITY'),但这只返回正确的'EPSG'。 EPSG是权威