2016-10-04 77 views
1

我正在尝试从gml文件中获取投影。这是文件的顶部:如何从gml文件中获取投影

<?xml version="1.0" encoding="UTF-8"?> <eop:Mask xmlns:eop="http://www.opengis.net/eop/2.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:gml="http://www.opengis.net/gml/3.2" gml:id="S2A_OPER_MSK_CLOUDS_SGS__20160914T145755_A006426_T31UCT_B00_MSIL1C"> <gml:name>MSK_CLOUDS pixels mask from data-strip S2A_OPER_MSK_CLOUDS_SGS__20160914T145755_A006426_T31UCT_B00_MSIL1C</gml:name> <gml:boundedBy> 
    <gml:Envelope srsName="urn:ogc:def:crs:EPSG:8.7:32631"> 
     <gml:lowerCorner>300000 5690220</gml:lowerCorner> 
     <gml:upperCorner>368340 5777580</gml:upperCorner> 
    </gml:Envelope> </gml:boundedBy> <eop:maskMembers> 
    <eop:MaskFeature gml:id="OPAQUE.0"> 
     <eop:maskType codeSpace="urn:gs2:S2PDGS:maskType">OPAQUE</eop:maskType> 
     <eop:extentOf> 
     <gml:Polygon gml:id="OPAQUE.0_Polygon"> 
      <gml:exterior> 
      <gml:LinearRing> 
       <gml:posList srsDimension="2">320340 5776020 320520 5776020 320520 5775960 320700 5775960 320700 5775900 320760 5775900 320760 5775840 320820 5775840 320820 5775660 320760 5775660 320760 5775600 320700 5775600 320700 5775540 320340 5775540 320340 5775600 320280 5775600 320280 5775660 320220 5775660 320220 5775900 320280 5775900 320280 5775960 320340 5775960 320340 5776020</gml:posList> 
      </gml:LinearRing> 
      </gml:exterior> 
     </gml:Polygon> 
     </eop:extentOf> 
    </eop:MaskFeature> 
... 

我试图使用从https://pcjericks.github.io/py-gdalogr-cookbook/projection.html代码:

from osgeo import ogr, osr 
driver = ogr.GetDriverByName('ESRI Shapefile') 
dataset = driver.Open(r'c:\data\yourshpfile.shp') 

# from Layer 
layer = dataset.GetLayer() 
spatialRef = layer.GetSpatialRef() 
# from Geometry 
feature = layer.GetNextFeature() 
geom = feature.GetGeometryRef() 
spatialRef = geom.GetSpatialReference() 

但spatialRef的两个版本都没有。

您可以从文件中看到投影是在边界框中给出的(它在第一行的末尾指出,然后是第2行中带有EPSG代码的信封)。(它没有'在文件的其他地方说'crs'或'EPSG')。

谁能告诉我如何访问投影信息?

我可能会以某种方式进入边界框然后获得投影吗?

回答

0

通常您可以在Granules提供的JPEG 2000文件中找到投影信息。

使用GDAL是:

gdalinfo是* .jp2

为您提供:

PROJCS["WGS 84/UTM zone 21N", 
GEOGCS["WGS 84", 
    DATUM["WGS_1984", 
     SPHEROID["WGS 84",6378137,298.257223563, 
      AUTHORITY["EPSG","7030"]], 
     AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich",0, 
     AUTHORITY["EPSG","8901"]], 
    UNIT["degree",0.0174532925199433, 
     AUTHORITY["EPSG","9122"]], 
    AXIS["Latitude",NORTH], 
    AXIS["Longitude",EAST], 
    AUTHORITY["EPSG","4326"]], 
PROJECTION["Transverse_Mercator"], 
PARAMETER["latitude_of_origin",0], 
PARAMETER["central_meridian",-57], 
PARAMETER["scale_factor",0.9996], 
PARAMETER["false_easting",500000], 
PARAMETER["false_northing",0], 
.... 

这就是所有我