Python遥感开发之GDAL读写遥感影像

• 1 读取tif信息方法一
• 2 读取tif信息方法二
• 3 自己封装读取tif的方法（推荐）
• 4 对读取的tif数据进行简单运算
• 5 写出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  # 查看波段个数，单波段默认是1# print(num_bands)# data_band = dataset.GetRasterBand(1)  # 1波段的具体内容# print(data_band.ReadAsArray())data = dataset.ReadAsArray()  # 转为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()#转为numpy格式data = data.astype(np.float32)#转为float类型a = data[0][0]data[data == a] = np.nan #原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值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数据进行简单运算

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()#转为numpy格式data = data.astype(np.float32)#转为float类型a = data[0][0]data[data == a] = np.nan #原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值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基础上所有值乘以2print(data[1])#可以进行条件筛选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()#转为numpy格式data = data.astype(np.float32)#转为float类型a = data[0][0]data[data == a] = np.nan #原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值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类型进行存储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基础上所有值乘以2#生成新的tifsave_tif(data,"lucc.tif","new_lucc.tif")#可以自己指定文件目录