本文记录了使用GDAL库读取遥感数据的代码，封装成了函数方便读写，只需要输入影像文件路径即可。并记录了使用numpy的计算过程的示例代码。
文章目录
- 一、环境配置
- 二、函数详解
- 三、读影像
- 四、写影像
- 五、应用示例
一、环境配置
import os
import numpy as np
from osgeo import gdal二、函数详解
1，GetGeoTransform()
GetGeoTransform()
GeoTransform[0]，左上角横坐标（应该是投影坐标）
GeoTransform[2]，行旋转
GeoTransform[1]，像元宽度(影像在水平空间的分辨率)
GeoTransform[3]，左上角纵坐标（应该是投影坐标）
GeoTransform[4]，列旋转
GeoTransform[5]，像元高度(影像在垂直空间的分辨率)，
如果影像是指北的，GeoTransform[2]和GeoTransform[4]这两个参数的值为0，GeoTransform[5]为负；
如果图像不含地理坐标信息，默认返回值是：(0,1,0,0,0,1)计算图像地理坐标
若图像中某一点的行数和列数分别为：row, column
则该点的地理坐标为：
xCoordinate = geoTransform[0] + column * geoTransform[1] + row * geoTransform[2]
yCoordinate = geoTransform[3] + column * geoTransform[4] + row * geoTransform[5]
2，ReadAsArray(xoff,yoff,xsize,ysize)
ReadAsArray(xoff,yoff,xsize,ysize)
#xoff,yoff是取值窗口的左上角在实际数据中所处象元的xy位置。
#xsize,ysize是取值窗口覆盖的区域大小
三、读影像
#读图像文件
def read_img(self,filename):
dataset=gdal.Open(filename) #使用gdal打开文件
im_width = dataset.RasterXSize #栅格矩阵的列数(每个波段的列数和)
im_height = dataset.RasterYSize #栅格矩阵的行数（每个波段的行数和）
im_geotrans = dataset.GetGeoTransform() #仿射矩阵，地理坐标信息
im_proj = dataset.GetProjection() #地图投影信息
im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵
del dataset
return im_proj,im_geotrans,im_data四、写影像
# 写文件，写成tiff
def write_img(filename, im_proj, im_geotrans, im_data):
# 判断栅格数据的数据类型
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 判读数组维数
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff") # 数据类型必须有，因为要计算需要多大内存空间
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
dataset.SetProjection(im_proj) # 写入投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset五、应用示例
import os
import numpy as np
from osgeo import gdal
class GRID:
#读图像文件
def read_img(self,filename):
dataset=gdal.Open(filename) #打开文件
im_width = dataset.RasterXSize #栅格矩阵的列数(每个波段的列数和)
im_height = dataset.RasterYSize #栅格矩阵的行数（每个波段的行数和）
im_geotrans = dataset.GetGeoTransform() #仿射矩阵
im_proj = dataset.GetProjection() #地图投影信息
im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵
del dataset
return im_proj,im_geotrans,im_data
#写文件，以写成tif为例
def write_img(self,filename,im_proj,im_geotrans,im_data):
#gdal数据类型包括
#gdal.GDT_Byte,
#gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
#gdal.GDT_Float32, gdal.GDT_Float64
#判断栅格数据的数据类型
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
#判读数组维数
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape
#创建文件
driver = gdal.GetDriverByName("GTiff") #数据类型必须有，因为要计算需要多大内存空间
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
dataset.SetProjection(im_proj) #写入投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) #写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
if __name__ == "__main__":
os.chdir(r'C:\Users\XX\Desktop\fx') #切换路径到待处理图像所在文件夹
run = GRID()
proj,geotrans,data = run.read_img('DEM.tif') #读DEM数据
print (proj)
print (geotrans)
print (data)
print (data.shape)
print (data.shape[0])
print (data.shape[1])
proj,geotrans,ndvi = run.read_img('NDVI.tif') #读NDVI数据
print (proj)
print (geotrans)
print (ndvi) #.astype(np.float)转换数值类型
print (ndvi.shape)
print (ndvi.shape[0])
print (ndvi.shape[1])
proj,geotrans,slope = run.read_img('slope.tif') #读slope数据
print (proj)
print (geotrans)
print (slope)
print (slope.shape)
print (slope.shape[0])
print (slope.shape[1])
#处理影像
divide = np.full((2226,2226),np.nan)
print(range(data.shape[1]))
for i in range(1,2225):
for j in range(1,2225):
if slope[i][j]>= 0.04366094 and ndvi[i][j]!=np.nan and data[i][j]!=np.nan and ndvi[i][j]!=0 and data[i][j]!=0:
divide[i][j]=ndvi[i][j]/data[i][j]
else:
divide[i][j]=np.nan
# ndvioutmin, ndvioutmax = ndviout.min(), ndviout.max() # 求最大最小值
# ndviout = (ndviout-ndvioutmin)/(ndvioutmax-ndvioutmin) # 归一化：(矩阵元素-最小值)/(最大值-最小值)
# print('=========elevation=========')
# print(divide.astype(np.double))
outResult = divide.astype(np.double)/1000
#导出影像
run.write_img('result.tif', proj, geotrans, outResult)
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
