遥感影像重投影

一、简介

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

二、重投影方法

　　GDAL是处理遥感影像强大的软件包，对遥感图像重投影，下面列举两种方法，

1 gdal.warp()

from osgeo import gdal
 
in_ds = gdal.Open(r'F:\algorithm\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

2 gdal.AutoCreateWarpedVRT()

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

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

重采样方法：

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：

from osgeo import gdal, osr
 
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
os.chdir(r'F:\algorithm')
old_ds = gdal.Open('lansat8_test.tif')
vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(),gdal.GRA_Bilinear)
gdal.GetDriverByName('gtiff').CreateCopy('landsat8_test_wgs84.tif', vrt_ds)

三、仿射变换

 

1 GetGeoTransform

#GetGeoTransform，在真实坐标和栅格数据坐标具有相同srs情况下，计算坐标偏移。
#作用：图像坐标（行列号）和现实世界坐标（投影坐标或地理坐标）之间的转换。是仿射变换，不是投影转换，和上面不同。
#0、3 x\y坐标 起始点现实世界坐标  1、5 像素宽度和高度  2、4 x\y方向旋转角
gt = ds.GetGeoTransform() #正变换：图像坐标到现实世界坐标。正变换时输入行列号，输出的现实世界坐标是栅格图像srs下的坐标.
inv_gt = gdal.InvGeoTransform(gt) #逆变换：现实世界坐标到图像坐标
offsets = gdal.ApplyGeoTransform(inv_gt, 465200, 5296000) #逆变换：输入的投影坐标具有和栅格图像的相同的srs
xoff, yoff = map(int, offsets)  #取整

2 gdal.Transformer

#gdal.Transformer，可计算相同srs下的坐标偏移；不能用于不同srs投影转换
#　　作用：现实世界坐标（投影坐标或地理坐标）与图像坐标（行列号）之间的转换、两个栅格图像之间像素坐标偏移（行列号），
如镶嵌#原理就是在相同srs情况下，计算图1的像素坐标到现实世界坐标的偏移，再从现实世界坐标偏移到图2的像素坐标。
其实就是两次仿射变换（正、逆），从而把图1的像素坐标偏移到图2的像素坐标。#所以不能用于不同srs情况，
因为该函数没有内置不同srs的投影转换公式。只能用于相同srs下，两个栅格数据集坐标的偏移。
#这里in_ds和out_ds具有相同srs。转换目的是为了把不同栅格数据的图像坐标（行列号）进行偏移，方便镶嵌
#in_ds是in_ds = gdal.Open()。
trans = gdal.Transformer(in_ds, out_ds, [])  #空白用于设置转换器选项
success, xyz = trans.TransformPoint(False, 0, 0) #False基于源数据计算目标栅格，反之为True。起始坐标为左上角 0，0
x,y,z = map(int, xyz)

#图像坐标和现实世界坐标之间转换
trans = gdal.Transformer(out_ds, None, [])
success, xyz = trans.TransformPoint(0, 1078, 648)

 

标签：投影,srs,遥感,影像,GRA,ds,gdal,坐标
来源： https://www.cnblogs.com/suoyike1001/p/15244791.html