0

我正在尝试通过 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 information from the shapefile to reproject the images to
outSpatialRef = lyr.GetSpatialRef().ExportToWkt()
ds = None # close file
print("\nSpatial referencing information of the shapefile:\n", outSpatialRef)

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

Spatial referencing information 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[2])
  box.AddPoint(extent[1], extent[3])
  box.AddPoint(extent[0], 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, 406749.40699999966, 565178.6085999999, 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))
4

1 回答 1

0

原来 shapefile 所在的坐标系非常关键，它应该在GCS_WGS_1984中。

于 2021-03-29T00:13:21.507 回答