Python遙感開發(fā)之GDAL讀寫遙感影像

• 1 讀取tif信息方法一
• 2 讀取tif信息方法二
• 3 自己封裝讀取tif的方法（推薦）
• 4 對讀取的tif數(shù)據(jù)進(jìn)行簡單運(yùn)算
• 5 寫出tif影像(推薦)

前言：主要介紹了使用GDAL讀寫遙感影像數(shù)據(jù)的操作，包括讀取行、列、投影、值以及數(shù)據(jù)的簡單運(yùn)算和生成新的tif影像。

---

1 讀取tif信息方法一

from osgeo import gdal import numpy as npif __name__ == '__main__':dataset = gdal.Open("lucc.tif")#讀取的是某地的土地類型col = dataset.RasterXSize # 圖像長度print("col:",col)row = dataset.RasterYSize # 圖像寬度print("row:", row)geotrans = dataset.GetGeoTransform() # 讀取仿射變換print("geotrans:", geotrans)proj = dataset.GetProjection() # 讀取投影print("proj:", proj)# num_bands = dataset.RasterCount # 查看波段個數(shù)，單波段默認(rèn)是1# print(num_bands)# data_band = dataset.GetRasterBand(1) # 1波段的具體內(nèi)容# print(data_band.ReadAsArray())data = dataset.ReadAsArray() # 轉(zhuǎn)為numpy格式data = data.astype(np.float32)a = data[0][0]data[data == a] = np.nanprint("data:", data)#遍歷每一行像元值for i in range(0,row):print(i,data[i])#遍歷讀取每一個像元for i in range(0,row):for j in range(0,col):if not np.isnan(data[i][j]):#篩選有效值print(data[i][j])

2 讀取tif信息方法二

from osgeo import gdalnumeric import numpy as npif __name__ == '__main__':data = gdalnumeric.LoadFile("lucc.tif")data = data.astype(np.float32)a = data[0][0]data[data == a] = np.nan#遍歷每一行像元for d in data:print(d)#遍歷每一個像元for d in data:for s in d:print(s)

3 自己封裝讀取tif的方法（推薦）

import numpy as np from osgeo import gdal,gdalnumericdef read_tif01(filepath):dataset = gdal.Open(filepath)col = dataset.RasterXSize#圖像長度row = dataset.RasterYSize#圖像寬度geotrans = dataset.GetGeoTransform()#讀取仿射變換proj = dataset.GetProjection()#讀取投影data = dataset.ReadAsArray()#轉(zhuǎn)為numpy格式data = data.astype(np.float32)#轉(zhuǎn)為float類型a = data[0][0]data[data == a] = np.nan #原因：讀取某一個行政區(qū)的影像圖的時(shí)候，往往第一行的第一列值為空值return [col, row, geotrans, proj, data]def read_tif02(filepath):data = gdalnumeric.LoadFile(filepath)data = data.astype(np.float32)a = data[0][0]data[data == a] = np.nanreturn dataif __name__ == '__main__':col, row, geotrans, proj, data = read_tif01("lucc.tif")data2 = read_tif02("lucc.tif")

4 對讀取的tif數(shù)據(jù)進(jìn)行簡單運(yùn)算

import numpy as np from osgeo import gdaldef read_tif01(filepath):dataset = gdal.Open(filepath)col = dataset.RasterXSize#圖像長度row = dataset.RasterYSize#圖像寬度geotrans = dataset.GetGeoTransform()#讀取仿射變換proj = dataset.GetProjection()#讀取投影data = dataset.ReadAsArray()#轉(zhuǎn)為numpy格式data = data.astype(np.float32)#轉(zhuǎn)為float類型a = data[0][0]data[data == a] = np.nan #原因：讀取某一個行政區(qū)的影像圖的時(shí)候，往往第一行的第一列值為空值return [col, row, geotrans, proj, data]if __name__ == '__main__':col, row, geotrans, proj, data = read_tif01("lucc.tif")print(data[1])data = data*2 #在原來的data基礎(chǔ)上所有值乘以2print(data[1])#可以進(jìn)行條件篩選data[data==6] = 10#所有像元值為6的重新賦值為10print(data[1])

5 寫出tif影像(推薦)

import numpy as np from osgeo import gdaldef read_tif01(filepath):dataset = gdal.Open(filepath)col = dataset.RasterXSize#圖像長度row = dataset.RasterYSize#圖像寬度geotrans = dataset.GetGeoTransform()#讀取仿射變換proj = dataset.GetProjection()#讀取投影data = dataset.ReadAsArray()#轉(zhuǎn)為numpy格式data = data.astype(np.float32)#轉(zhuǎn)為float類型a = data[0][0]data[data == a] = np.nan #原因：讀取某一個行政區(qū)的影像圖的時(shí)候，往往第一行的第一列值為空值return [col, row, geotrans, proj, data]def save_tif(data, file, output):ds = gdal.Open(file)shape = data.shapedriver = gdal.GetDriverByName("GTiff")dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Float32)#以float類型進(jìn)行存儲dataset.SetGeoTransform(ds.GetGeoTransform())dataset.SetProjection(ds.GetProjection())dataset.GetRasterBand(1).WriteArray(data)if __name__ == '__main__':col, row, geotrans, proj, data = read_tif01("lucc.tif")data = data*2 #在原來的data基礎(chǔ)上所有值乘以2#生成新的tifsave_tif(data,"lucc.tif","new_lucc.tif")#可以自己指定文件目錄

總結(jié)

以上是生活随笔為你收集整理的Python遥感开发之GDAL读写遥感影像的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。