Python地学分析 — GDAL对遥感影像重投影

欢迎关注博主的微信公众号：“智能遥感”。

该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。

本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

“ 对遥感影像重投影是遥感数据预处理的常见手段之一，本节通过gdalwarp和建立VRT分别对遥感影像重投影。”

今天的遥感之美封面图 — 八城外有回城处，哈密伊犁吐鲁番。一提到吐鲁番，是不是想到王翰笔下的葡萄美酒夜光杯，抑或是《西游记》中“喋喋不休”的火焰山。现在让我们通过太空之眼，来尽情地俯瞰这如梦似幻的美景吧。

该图像由Landsat 7的Enhanced Thematic Mapper plus（ETM +）传感器获得（2017年8月7日）。

吐鲁番盆地位于我国新疆天山东部南坡，是一个奇特的盐湖与沙丘混合体。它是世界上少数几个位于海平面以下的景观之一。吐鲁番盆地是一个典型的地堑盆地，也是我国地势最低（-154.31米）和夏季气温最高的地方。大部分地面在海拔500米以下。

图像重投影

本节重点内容是对栅格数据进行重新投影，但它比矢量数据更复杂。对于矢量，你只需要每个顶点的新坐标就可以轻松实现投影转换，但对于栅格，你需要处理像元发生形变和偏移的情况，以及从旧单元格位置到新单元格位置的一对一映射。新单元格位置不存在（如下图）。确定新单元格像元值的最简单方法是使用最接近输出单元格映射的输入单元格中的值，这称为最近邻，是最快的方法，通常是分类数据所需的方法。除此之外的所有其他采样类型都将更改像元值，对于数据分类而言，这是我们不希望看到的结果。但是，如果使用最近邻方法，连续数据的栅格通常看起来不太美观，也就是图像上的物体可能出现“断裂”。对于这种情况，我们通常使用双线性插值或三次卷积，它采用的是周围像元的平均值。确定如何使用哪种插值方法，只能具体问题具体分析啦。

投影栅格时像元如何移动的示例。三角形和圆形是两个不同栅格的像元中心点。三角形是从圆圈来自的栅格的重新投影所创建的，注意图像尺寸甚至不同。

GDAL是处理遥感影像强大的软件包，对遥感图像重投影，下面列举两种方法，供大家参考，对每种方法均有测试，并进行了对比，小编发现，不同的方法得出的结果还是有差别的，例如同一个像元值会略有不同，或者整幅图像的行列数也会有所差异，各位小伙伴们做的时候需要注意一下，但这几种方法均可以使用。

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/11/26 22:40'# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/11/26 22:40'from osgeo import gdalin_ds = gdal.Open('D:/MOD09A1.A2017361.h28v06.006.2018005034659.hdf')
# 返回结果是一个list，list中的每个元素是一个tuple，每个tuple中包含了对数据集的路径，元数据等的描述信息
# tuple中的第一个元素描述的是数据子集的全路径
datasets = in_ds.GetSubDatasets()# 取出第1个数据子集（MODIS反射率产品的第一个波段）进行转换
# 第一个参数是输出数据，第二个参数是输入数据，后面可以跟多个可选项
gdal.Warp('D:/reprojection01.tif', datasets[0][0], dstSRS='EPSG:32649')  # UTM投影
gdal.Warp('D:/reprojection02.tif', datasets[0][0], dstSRS='EPSG:4326')   # 等经纬度投影# 关闭数据集
root_ds = None

下面是几种结果图：

UTM投影下的MOD09A1地表反射率图像

等经纬度投影下的MOD09A1地表反射率图像

NASA官网等经纬度影下的MOD09A1地表反射率图像

除了使用GDAL附带的gdalwarp实用程序之外，重新投影图像的最简单方法是使用VRT。有一个方便的功能，当你提供空间参考信息时，会为你创建重新投影的VRT数据集。具体参数如下：

src_ds是要重新投影的数据集。
src_wkt是源空间参考系的WKT表示。默认值为None，在这种情况下，它将使用源栅格中的SRS信息。如果此栅格没有SRS信息，则需要在此处提供。如果使用None，你也可以在这里提供。
dst_wkt是所需空间参照系的WKT表示。默认值为None，在这种情况下不会执行重新投影操作。
eRasampleAlg是表1中的重采样方法之一。默认值为GRA_NearestNeighbour。
maxerror是你要允许的最大错误量（以像元为单位）。默认值为0，用于精确计算。

表1 重采样方法

常量	描述
GRA_NearestNeighbour	选取最邻近的像元
GRA_Bilinear	邻近4个像元加权平均
GRA_Cubic	邻近的16个像元平均
GRA_CubicSpline	16个像元的三次B样条
GRA_Lanczos	36个像元Lanczos窗口
GRA_Average	求均值
GRA_Mode	出现频率最多的像元值

AutoCreateWarpedVRT函数不会在磁盘上创建VRT文件，但会返回一个数据集对象，然后可以使用CreateCopy将其保存为其他格式。以下示例采用使用UTM空间参考的自然颜色Landsat图像，创建具有等经纬度投影WGS84的目标SRS的扭曲VRT，并将VRT复制到GeoTIFF：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/11/26 22:40'# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/11/26 22:40'from osgeo import gdal, osrsrs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
os.chdir(r'D:\osgeopy-data\Landsat\Washington')
old_ds = gdal.Open('nat_color.tif')
vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(),gdal.GRA_Bilinear)
gdal.GetDriverByName('gtiff').CreateCopy('nat_color_wgs84.tif', vrt_ds)

该程序输出如下图所示。

使用UTM投影的原始Landsat图像，以及使用等经纬度投影的新图像

Python地学分析 — GDAL对遥感影像重投影相关推荐

Python地学分析 — GDAL将多个遥感图像叠加保存为tif文件
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
python/gdal处理遥感影像（读取、投影转换、裁剪、建立图像金字塔等）
python/gdal处理遥感影像(读取.投影转换.裁剪.建立图像金字塔等) gdal库简单介绍 python使用gdal 一.安装python环境二.安装gdal库三.使用gdal处理遥感影像 ...
gdal进行遥感影像读写_如何使用遥感影像进行矿物勘探
gdal进行遥感影像读写 Meet Jose Manuel Lattus, a geologist from Chile. In the latest Soar Cast, he discusses ...
Python地学分析 — 建立矢量数据缓冲区 06
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. Python的小伙伴们,你们好!上一 ...
Python地学分析 — 地理空间参考系介绍
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
Python地学分析 — 矢量数据集介绍
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 什么叫矢量数据,在这里,小编先给大家 ...
原创程序|基于GDAL的遥感影像批量处理工具介绍（三）
本文主要介绍基于C#+GDAL-Python实用工具开发的遥感影像批量处理工具,该工具目前主要包括影像批量切片生成KML文件和影像批量切片生成Tiles文件.该工具.Net框架版本为4.0,GDAL版 ...
8影像计算ndvi landsat_使用GDAL读取遥感影像的信息
读取影像数据集的元数据 GDAL已经提供了足够方便的函数,可以读取影像的一些元数据信息, 从而方便对数据进行处理.GDAL一般是以字典的形式对元数据进行组织的, 但是对于不同的栅格数据类型,元数据的类 ...
python整形怎么切片_遥感影像切分切片
遥感影像切片生活当中,我们可能经常遇到处理一个很大的遥感影像情况,例如做一些逐像元运算,不便于并行处理.因此制作了将影像切分成多个小片的程序,切分后保持原有的数值.数据类型.波段数.投影. 切片用 ...

Python地学分析 — GDAL对遥感影像重投影

Python地学分析 — GDAL对遥感影像重投影相关推荐

最新文章

热门文章