一、簡介

　　柵格資料進行重新投影比向量資料更復雜，對於向量，你只需要每個頂點的新座標就可以輕鬆實現投影轉換，但對於柵格，你需要處理像元發生形變和偏移的情況，以及從舊單元格位置到新單元格位置的一對一對映。新單元格位置不存在（如下圖）。確定新單元格像元值的最簡單方法是使用最接近輸出單元格對映的輸入單元格中的值，這稱為最近鄰，是最快的方法，通常是分類資料所需的方法。除此之外的所有其他取樣型別都將更改像元值，對於資料分類而言，這是我們不希望看到的結果。但是，如果使用最近鄰方法，連續資料的柵格通常看起來不太美觀，也就是影象上的物體可能出現“斷裂”。對於這種情況，我們通常使用雙線性插值或三次卷積，它採用的是周圍像元的平均值。

二、重投影方法

　　GDAL是處理遙感影像強大的軟體包，對遙感影象重投影，下面列舉兩種方法，

1gdal.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)

三、仿射變換

1GetGeoTransform

#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)  #取整

2gdal.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)