【Python+GIS】圖像轉置、定義栅格投影

圖像轉置、定義栅格投影

• 1. 實驗背景與資料
• 2. 實驗說明
• • 2.1 實驗目的
  • 2.2 參考文獻
• 3.代碼
• • 3.1 轉置
  • 3.2 定義投影
• 4. 總結與反思

1. 實驗背景與資料

最近下載下傳了程潔老師的LST資料，這個資料集填充了MOD11 LST的空缺，遂使用此資料集進行後面的操作。

資料可以在國家青藏高原資料中心開放擷取，貼一張資料簡介：

此資料集為hdf格式，包含有4個波段圖層，緯度、經度、LST和LST的品質，但是這樣分開就有一個小問題，把LST拖進ArcGIS裡面，會顯示沒有投影，而且并不是正的，如下圖所示：

也就是說，咱們想要使用這個資料，得進行一些小處理。

2. 實驗說明

2.1 實驗目的

（1）轉置。将豎着的圖像轉為橫着。

（2）定義投影。加入經緯度地理坐标。

2.2 參考文獻

《GDAL(python) 之GeoTransform》

《【GIS】GIS中仿射變換的六參數》

《python+GDAL給圖像設定地理參考（GeoTransform和Projection的設定）》

3.代碼

3.1 轉置

from osgeo import gdal, gdalconst
import numpy as np
import sys
import copy
import os

outputFolder=r"C:\Users\10366\Desktop\test\test\transfer\2"
inputFloder1=r"C:\Users\10366\Desktop\test\test\output"

input_folder_list_LN=os.listdir(inputFloder1) #讀取檔案夾裡所有檔案
LN_list=[] #建立一個隻裝tif格式的清單
for filename in input_folder_list_LN:  #周遊 
    filelast=filename[-3:]
    if  filelast=="tif":
        LN_list.append(filename) 

referencepath=r"C:\Users\10366\Desktop\test\LST REFERENCE.tif"

for i in range(0,len(LN_list)):
    img_path=inputFloder1+"\\"+LN_list[i]
    dataset = gdal.Open(img_path)  # 取出第j個波段
    img_width = dataset.RasterXSize
    img_height = dataset.RasterYSize
    img_data = np.array(dataset.ReadAsArray(0, 0, img_width, img_height), dtype=float)  # 将資料寫成數組，對應栅格矩陣
    img_data=img_data.T  #轉置！

# 儲存為TIF格式
    save_path=outputFolder+"\\"+LN_list[i]
    driver = gdal.GetDriverByName("GTiff")
    datasetnew = driver.Create(save_path, img_height,img_width, 1, gdal.GDT_Float32)
    #datasetnew.SetGeoTransform(adf_GeoTransform)
    #datasetnew.SetProjection(im_Proj)
    band = datasetnew.GetRasterBand(1)
    band.WriteArray(img_data)
    datasetnew.FlushCache()  # Write to disk.必須有清除緩存
    
print("OK")
           

轉置的重點是T，這在矩陣中是較為常見的：

3.2 定義投影

from osgeo import gdal, osr
from osgeo import gdalconst
import sys
import copy
import os
import numpy as np

inputFloder1=r"C:\Users\10366\Desktop\test\test\transfer\initial"
outputFolder=r"C:\Users\10366\Desktop\test\test\transfer\3"
reference=r"C:\Users\10366\Desktop\test\LST REFERENCE.tif" #參考影像


#LN檔案夾
input_folder_list_LN=os.listdir(inputFloder1) #讀取檔案夾裡所有檔案
LN_list=[] #建立一個隻裝tif格式的清單
for filename in input_folder_list_LN:  #周遊 
    filelast=filename[-3:]
    if  filelast=="tif":
        LN_list.append(filename)  

for i in range(0,len(LN_list)):
    LN=gdal.Open(inputFloder1+"\\"+LN_list[i])
    ref=gdal.Open(reference)  

    #geotransform = LN.GetGeoTransform()
    projection=ref.GetProjection() #參考影像的投影坐标系
    cols = LN.RasterXSize
    rows = LN.RasterYSize
    
    LNBand = LN.GetRasterBand(1)
    LNData = LNBand.ReadAsArray(0,0,cols,rows)
    
    xmin=73.55  #經度最小值
    xmax=134.99 #經度最大值
    ymin=18.33 #緯度最小值
    ymax=53.47 #緯度最大值
    xscale=(xmax-xmin)/cols #經度方向的分辨率
    yscale=(ymax-ymin)/rows #緯度方向的分辨率
    GeoList=[xmin,xscale,0.0,ymax,0.0,-yscale] #設定地理坐标轉換參數
    geotransform=tuple(GeoList) #建立地理坐标轉換

    result = LNData 
    resultPath = outputFolder+"\\pro1.tif" #結果的儲存路徑名
    print(resultPath)

    format = "GTiff"   
    driver = gdal.GetDriverByName(format)
    ds = driver.Create(resultPath, cols, rows, 1, gdal.GDT_Float32)
    ds.SetGeoTransform(geotransform)
    ds.SetProjection(projection)
    band  = ds.GetRasterBand(1)
    band.WriteArray(result)
    ds.FlushCache()  
    ds = None

print("ok")
           

結果如下圖所示：

4. 總結與反思

1. 轉置時，由于數組的shape變了，是以在建立影像的時候，要注意行數和列數的設定。

img_width = dataset.RasterXSize #原來的列數
img_height = dataset.RasterYSize #原來的行數
datasetnew = driver.Create(save_path, img_height,img_width, 1, gdal.GDT_Float32) 
           

1. geotransform是元組，元組的資料是不能修改的，是以得先建立list，再list轉元組tuple

GeoList=[xmin,xscale,0.0,ymax,0.0,-yscale] #設定地理坐标轉換參數
geotransform=tuple(GeoList) #建立地理坐标轉換
           

1. 在設定地理坐标轉換參數的時候，有六個參數，分别是：

   （左上角x坐标， 水準分辨率，旋轉參數， 左上角y坐标，旋轉參數，豎直分辨率）

   **要注意前面是先分辨率再旋轉參數，後面是先旋轉參數再分辨率，**我之前就是以為這倆是一樣的，屢戰屢敗，圖像根本出不來。

完成上述操作後，就處理好了，可以使用了。