目前从Google、Bing等GIS平台可获取免费全球高分辨率卫星影像，但其官方软件（如Google Earth）通常只支持单幅影像的保存，无法批量获取一个区域内的所有影像。国内也有91卫图、水经注等下载助手，通过转发可获得免费权限，可以批量获取低分辨率影像，而高分辨率影像则收费。博主在使用QGIS时，发现了一个可以批量获取Google、Bing等各级分辨率卫星影像的方法。
该方法的基本原理如下：QGIS是一个功能丰富的开源GIS软件，它可以方便地导入Google、Bing等卫星影像底图，并提供Python的API接口，可以通过代码将指定坐标范围的底图保存为图片和空间参考文件，写几个循环就可以实现批量获取了。
目录
一、QGIS安装
QGIS是一个类似于ArcGIS的开源GIS软件，个人使用体验非常不错，免费无需破解，多平台（Windows、Linux、macOS等），更新快，功能齐全，而且大数据的渲染比ArcGIS快很多，整体肯定还是不如ArcGIS（毕竟是商业软件），但完全可以满足使用，而且有很多非常有特色且实用的功能，这次讲到的方法就要用到，感兴趣的朋友可以试试！
其下载地址如下：Download QGIS。
有两种安装方式：1.（不推荐）使用OSGeo4W下载器（全家桶，里面还有Grass等其他软件可用）下载安装，比较麻烦，且由于是外网服务器，下载很容易失败；2.（推荐）直接下载安装包安装，提供最新版和稳定版，凭喜好都可以。接下来的演示以QGIS3.14.1为例。
二、导入卫星影像底图
打开QGIS，有两个方法导入底图：详见在QGIS中导入Google、Bing等地图和卫星影像。
这里以Bing卫星影像为例。
该类型的底图投影通常为EPSG:3857 - WGS 84 / Pseudo-Mercator。
底部状态栏，coordinate显示的是当前鼠标坐标，点击红色框中的按钮，可转换为当前视图的坐标范围。
对底图进行计算前要确定当前project的投影是否与底图一致，后面提供的代码进行的计算是以该投影为基础的，最好在打开QGIS后，首先打开底图，这样就是3857投影，如果是其他投影，坐标就得自己另外进行换算。
三、通过Python接口实现影像批量获取
打开Python Console，点击Show Editor，即可输入代码。
基本思路为：
1.选定分辨率（默认为96 dpi）和比例尺，并计算像素高宽（即单个像素在该空间坐标系下的大小）;
2.计算待裁剪影像的像素尺寸和坐标范围（确定的分辨率和比例尺下，像素尺寸与坐标范围是相互对应的，即要么根据指定的像素尺寸计算坐标范围，要么根据指定的坐标范围计算像素尺寸）；
3.设置输出参数，保存为图片。
PS：通过循环可以实现批量获取。
0.像素高宽的计算
首先说明下如何计算像素高宽，导入底图，随便选择一个区域，选择菜单Project→Import/Export→Export Map to Image…。
指定比例尺Scale（以1:10000为例）和分辨率Resolution（默认为96 dpi）。
根据Extent和Output Width/Height计算像素高宽。参考代码如下：
# -*- coding: utf-8 -*-
# 坐标范围[西, 南, 东, 北]
extent = [5730688.2950, 5299790.0539, 5734969.2534, 5302136.9080]
# 输出像素尺寸[高, 宽]
output_size = [886, 1618]
px_geosize = [(extent[3] - extent[1]) / output_size[0], (extent[2] - extent[0]) / output_size[1]]
print('Geosize of Pixel:', px_geosize)
输出结果为：
[2.645833374536395, 2.6488195259594955]
得到像素高宽后可以开始裁剪，这里提供两个不同裁剪方式的参考代码，有其他需求的可以在这基础上修改，需要用到Python的GDAL包，QGIS有内置，无需另外下载安装。
1.获取每个Polygon外接矩形范围的影像（固定比例尺和分辨率）
PS：注意Shapefile的投影与底图一致。
# -*- coding: utf-8 -*-
import os
from osgeo import ogr
shp_path = r'F:\GIS\Samples\Shapefile\Test.shp'
img_path = r'F:\GIS\Samples\PNG'
img_suffix = 'png'
# 像素[高, 宽]
px_geosize = [2.645859085290482, 2.6458015267176016]
# 在外接矩形的[上, 下, 左, 右]扩展的像素尺寸，无需扩展则为0
pad_size = [150, 150, 150, 150]
pad_geosize = [a * b for a, b in zip(pad_size, px_geosize)]
shp = ogr.Open(shp_path)
lyr = shp.GetLayer()
for feat in lyr:
# 默认Shapefile存在Name字段，若无请自行定义名称
name = str(feat.GetField('Name'))
geom = feat.GetGeometryRef()
geom_ext = geom.GetEnvelope()
# 待裁剪影像的坐标范围[min_x, min_y, max_x, max_y]
clip_ext = [geom_ext[0] - pad_geosize[0], geom_ext[2] - pad_geosize[1],
geom_ext[1] + pad_geosize[2], geom_ext[3] + pad_geosize[3]]
# 待裁剪影像的像素尺寸
clip_size = [int((clip_ext[2] - clip_ext[0]) / px_geosize[0]), int((clip_ext[3] - clip_ext[1]) / px_geosize[1])]
# 将待裁剪影像的坐标范围转为QGIS格式
rect = QgsRectangle(clip_ext[0], clip_ext[1], clip_ext[2], clip_ext[3])
# 图片保存设置
settings = iface.mapCanvas().mapSettings()
# 设置图片分辨率，默认为96 dpi，若不改可以不写
settings.setOutputDpi(96)
# 设置坐标范围
settings.setExtent(rect)
# 设置像素尺寸
settings.setOutputSize(QSize(clip_size[0], clip_size[1]))
job = QgsMapRendererSequentialJob(settings)
job.start()
job.waitForFinished()
image = job.renderedImage()
image.save(os.path.join(img_path, name + '.' + img_suffix))
print('task finished!')
2.获取Polygon范围内的影像（固定比例尺、分辨率和像素尺寸）
PS：注意Shapefile的投影与底图一致。
# -*- coding: utf-8 -*-
import os
from osgeo import ogr
shp_path = r'F:\GIS\Samples\Shapefile\Test.shp'
img_path = r'F:\GIS\Samples\PNG'
img_suffix = 'png'
# 像素[高, 宽]
px_geosize = [2.645859085290482, 2.6458015267176016]
# 单幅影像的像素尺寸
clip_size = [500, 500]
# 重叠区域的像素尺寸
overlay_size = [100, 100]
# 单幅影像的坐标尺寸
clip_geosize = [a * b for a, b in zip(clip_size, px_geosize)]
# 重叠区域的坐标尺寸
overlay_geosize = [a * b for a, b in zip(overlay_size, px_geosize)]
# 每次裁剪偏移的坐标尺寸
offset_geosize = [a - b for a, b in zip(clip_geosize, overlay_geosize)]
shp = ogr.Open(shp_path)
lyr = shp.GetLayer()
for feat in lyr:
# 默认Shapefile存在Name字段，若无请自行定义名称
polygon_name = str(feat.GetField('Name'))
geom = feat.GetGeometryRef()
# Polygon坐标范围[min_x, max_x, min_y, max_y]
geom_ext = geom.GetEnvelope()
# 以重叠尺寸略微扩大坐标范围，可删除
geom_ext = [geom_ext[0] - overlay_geosize[1], geom_ext[1] + overlay_geosize[1],
geom_ext[2] - overlay_geosize[0], geom_ext[2] + overlay_geosize[0]]
# 坐标范围内的按裁剪尺寸划分的条带数量
clip_count = [int((geom_ext[3] - geom_ext[2]) / offset_geosize[0]) + 1,
int((geom_ext[1] - geom_ext[0]) / offset_geosize[1]) + 1]
for i in range(clip_count[0]):
for j in range(clip_count[1]):
# 待裁剪影像的坐标范围[min_x, min_y, max_x, max_y]
clip_ext = (geom_ext[0] + offset_geosize[0] * j,
geom_ext[3] - offset_geosize[0] * i - clip_geosize[0],
geom_ext[0] + offset_geosize[0] * j + clip_geosize[1],
geom_ext[3] - offset_geosize[0] * i)
# 根据单幅影像的地理范围创建Polygon
clip_ring = ogr.Geometry(ogr.wkbLinearRing)
clip_ring.AddPoint(clip_ext[0], clip_ext[1])
clip_ring.AddPoint(clip_ext[2], clip_ext[1])
clip_ring.AddPoint(clip_ext[2], clip_ext[3])
clip_ring.AddPoint(clip_ext[0], clip_ext[3])
clip_ring.AddPoint(clip_ext[0], clip_ext[1])
clip_geom = ogr.Geometry(ogr.wkbPolygon)
clip_geom.AddGeometry(clip_ring)
# 若该幅影像坐标范围不与任何Polygon相交，则跳过裁剪
if not geom.Intersect(clip_geom):
continue
# 设置名称
name = polygon_name + '_{0}_{1}'.format(str(i).zfill(3), str(j).zfill(3))
# 将待裁剪影像的坐标范围转为QGIS格式
rect = QgsRectangle(clip_ext[0], clip_ext[1], clip_ext[2], clip_ext[3])
# 图片保存设置
settings = iface.mapCanvas().mapSettings()
# 设置图片分辨率，默认为96 dpi，若不改可以不写
settings.setOutputDpi(96)
# 设置坐标范围
settings.setExtent(rect)
# 设置像素尺寸
settings.setOutputSize(QSize(clip_size[1], clip_size[0]))
job = QgsMapRendererSequentialJob(settings)
job.start()
job.waitForFinished()
image = job.renderedImage()
image.save(os.path.join(img_path, name + '.' + img_suffix))
print('task finished!')
四、空间参考文件（World File）的生成
World File为Esri公司制定的给普通栅格图片添加空间参考的纯文本文件格式，可以使用记事本直接打开。World File文件名和主图片文件名一致，后缀名一般为3个字符，由主图片后缀名缩写加上字符w（如pgw、tfw、jgw等）。虽后缀名各有不同，但其读写方式与其他纯文本文件格式（如txt）一致。参考代码如下：
# -*- coding: utf-8 -*-
import os
# 像素[高, 宽]
px_geosize = [2.645859085290482, 2.6458015267176016]
# 影像的坐标范围[min_x, min_y, max_x, max_y]
clip_ext = [-7966406.0, 1203525.0, -7962125.1, 1205871.9]
with open(os.path.join(wf_path, name + '.' + wf_suffix), 'w') as f:
f.writelines([str(px_geosize[0]) + '\n', str(0) + '\n', str(0) + '\n',
str(-px_geosize[1]) + '\n', str(clip_ext[0]) + '\n', str(clip_ext[3]) + '\n'])
常用地学数据的空间参考相关信息可参考: GeoTIFF、Shapefile和World File中空间参考的获取及其参数意义。
912
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
