有没有办法在 Python 中将多边形 shapefile 转换为坐标？

问题描述

我正在尝试通过 ESA Sentinel 数据中心从 Sentinel 2 下载卫星图像。

我用来获取 shapefile 图层范围以设置查询的代码不是经纬度坐标，而是奇怪的数字。我仔细地遵循了实际说明，但没有运气。

任何有关如何解决此问题的建议或帮助将不胜感激！

代码如下：

# Get the shapefile layer's extent
driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(shapefile,0)
lyr = ds.GetLayer()
extent = lyr.GetExtent()
print("Extent of the area of interest (shapefile):\n",extent)

# get projection @R_914_4045@ion from the shapefile to reproject the images to
outSpatialRef = lyr.GetSpatialRef().ExportToWkt()
ds = None # close file
print("\nSpatial referencing @R_914_4045@ion of the shapefile:\n",outSpatialRef)

Extent of the area of interest (shapefile):
 (363337.9978,406749.40699999966,565178.6085999999,633117.0013999995)

Spatial referencing @R_914_4045@ion of the shapefile:
 PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB 1936",DATUM["OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],ParaMETER["false_easting",400000.0],ParaMETER["false_northing",-100000.0],ParaMETER["central_meridian",-2.0],ParaMETER["scale_factor",0.9996012717],ParaMETER["latitude_of_origin",49.0],UNIT["Meter",1.0]]

#   extent of our shapefile in the right format for the Data Hub API.
def bBox(extent):
  # Create a polygon from the extent tuple
  Box = ogr.Geometry(ogr.wkbLinearRing)
  Box.AddPoint(extent[0],extent[2])
  Box.AddPoint(extent[1],extent[3])
  Box.AddPoint(extent[0],extent[2])
  poly = ogr.Geometry(ogr.wkbpolygon)
  poly.AddGeometry(Box)
  return poly

# Let's see what it does
print(extent)
print(bBox(extent))

(363337.9978,633117.0013999995)
polyGON ((363337.9978 565178.6086 0,406749.407 565178.6086 0,406749.407 633117.001399999 0,363337.9978 633117.001399999 0,363337.9978 565178.6086 0))

解决方法

事实证明 shapefile 所在的坐标系非常重要，它应该在 GCS_WGS_1984 中。