用shape多边形文件下载对应的谷歌地球高清影像
前言
我们想要利用遥感图像进行处理分析,但是单张下载起来很麻烦,我们可以根据对应的shape多边形文件,实现批量化下载。
一、准备好shape文件数据
二、在QGIS的Python编辑器中复制以下代码
代码如下(示例)2种方法:
1、下载MultiPolygon对象中每个子MultiPolygon的外接矩形
import processing
import geopandas as gpd
import numpy as np
import os
from shapely.geometry import Polygon
# 选择目标shp文件,注意投影坐标系要和底图一致
data = gpd.read_file('G:/桌面/data append/China_PV_training_polygon_2020.shp')
image_paved_out_dir = 'G:/桌面/新建文件夹 (4)/'
for i, geom in enumerate(data['geometry']):
# 检查是否为多边形或多边形集合
if isinstance(geom, Polygon):
polygons = [geom]
elif geom.is_empty:
continue # 如果几何对象为空,则跳过
else:
polygons = list(geom.geoms) # 从MultiPolygon获取所有Polygon
for poly in polygons:
xx, yy = poly.exterior.coords.xy
extent_0 = np.min(xx)
extent_1 = np.max(xx)
extent_2 = np.min(yy)
extent_3 = np.max(yy)
extent = f"{extent_0},{extent_1},{extent_2},{extent_3} [EPSG:4326]"
# 测试:打印每个多边形的地理范围
print(f"多边形 {i}: 范围 = {extent}")
out_path = os.path.join(image_paved_out_dir, f"{i}.tif")
# 下载影像
processing.run("native:rasterize", {
'EXTENT': extent,
'EXTENT_BUFFER': 0,
'TILE_SIZE': 512,
'MAP_UNITS_PER_PIXEL': 0.000005,
'MAKE_BACKGROUND_TRANSPARENT': False,
'MAP_THEME': None,
'LAYERS': ['type=xyz&zmin=0&zmax=20&url=https://mt1.google.com/vt/lyrs%3Ds%26x%3D{x}%26y%3D{y}%26z%3D{z}'],
'OUTPUT': out_path
})
2、下载一个MultiPolygon对象的外接矩形
import processing
import geopandas as gpd
import numpy as np
import os
from shapely.geometry import Polygon, MultiPolygon
# 选择目标shp文件
data = gpd.read_file('G:/桌面/data append/China_PV_training_polygon_2020.shp')
image_paved_out_dir = 'G:/桌面/新建文件夹 (4)/'
for i, geom in enumerate(data['geometry']):
# 创建空列表以存储所有多边形的坐标
all_x_coords = []
all_y_coords = []
# 检查几何对象类型并相应处理
if isinstance(geom, (Polygon, MultiPolygon)):
if isinstance(geom, Polygon):
polygons = [geom]
else:
polygons = list(geom.geoms)
# 提取所有子多边形的坐标
for poly in polygons:
xx, yy = poly.exterior.coords.xy
all_x_coords.extend(xx)
all_y_coords.extend(yy)
# 计算整个MultiPolygon的范围
extent_0 = np.min(all_x_coords)
extent_1 = np.max(all_x_coords)
extent_2 = np.min(all_y_coords)
extent_3 = np.max(all_y_coords)
extent = f"{extent_0},{extent_1},{extent_2},{extent_3} [EPSG:4326]"
out_path = os.path.join(image_paved_out_dir, f"{i}.tif")
# 下载影像
processing.run("native:rasterize", {
'EXTENT': extent,
'EXTENT_BUFFER': 0,
'TILE_SIZE': 512,
'MAP_UNITS_PER_PIXEL': 0.000005,
'MAKE_BACKGROUND_TRANSPARENT': False,
'MAP_THEME': None,
'LAYERS': ['type=xyz&zmin=0&zmax=20&url=https://mt1.google.com/vt/lyrs%3Ds%26x%3D{x}%26y%3D{y}%26z%3D{z}'],
'OUTPUT': out_path
})
elif geom.is_empty:
continue # 跳过空的几何对象