基于python的gdal读取遥感影像(基于Python的图书馆信息管理系统研发)

基于python的gdal读取遥感影像 基于python的gdal读写遥感影像1. gdal介绍2. 代码详解2.1 读取数据2.2 写入影像3. 完整案例1. gdal介绍

推荐整理分享基于python的gdal读取遥感影像(基于Python的图书馆信息管理系统研发)，希望有所帮助，仅作参考，欢迎阅读内容。

文章相关热门搜索词:基于Python的逻辑回归算法,基于Python的招聘网站数据分析,基于Python的电影推荐系统,基于Python的大数据分析,基于Python的旅游推荐系统,基于Python的金融分析与风险管理(第2版)pdf,基于python的GDP增长模型,基于Python的大数据分析,内容如对您有帮助，希望把文章链接给更多的朋友！

GDAL(Geospatial Data Abstraction Library)主要用来读取地理空间数据，现在的GDAL包并不是单独的GDAL，而是集成了GDAL和OGR的。OGR用于处理矢量数据。因此，GDAL既可以用来处理栅格也可以处理矢量文件。

2. 代码详解

下面代码主要为读写tif的代码，图像其实就是一个二维矩阵，其他格式的遥感影像都可以参考。包括nc，hdf等等。

2.1 读取数据# 调包from osgeo import gdalimport numpy as np打开影像：gdal.Open()为打开图像的操作，gdal.GetDriverByName()为注册哪种格式的数据。通常影像为Geotiff格式，也可以不写参数，默认注册所有格式。# 打开图像并创建空间 data = gdal.Open(filename) driver = gdal.GetDriverByName('GTiff')获取数据的基本信息：# 数据集的基本信息 print('Raster Driver : {d}\n'.format(d=driver.ShortName)) print('影像的波段数: ', data.RasterCount) img_width, img_height = data.RasterXSize, data.RasterYSize print('影像的列，行数: {r}rows * {c}colums'.format(r=img_width, c=img_height)) print('栅格数据的空间参考：{}'.format(data.GetGeoTransform())) # 栅格数据的6参数 print('投影信息：{}\n'.format(data.GetProjection())) # 栅格数据的投影

输出如下： 关于数据的基本信息，其中数据的空间参考和投影信息最为重要。在写影像时是必不可少的参数。GetGeoTransform()为获取数据的仿射变换矩阵，6个参数，记录了影像的左上角坐标，旋转，xy方向上的分辨率；GetProjection()为获取影响的投影信息，即影像的地理参考，一般都为WGS84。 若数据格式为NC或者HDF，仿射变换矩阵需要计算。

Raster Driver : GTiff影像的波段数: 6影像的列行数: 1143rows * 721colums栅格数据的空间参考：(429255.0, 30.0, 0.0, 4429725.0, 0.0, -30.0)投影信息：PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]将数据转换为数组形式，方便进行后期处理。 注意：python中的索引从0开始到n-1，而波段读取的索引从1到n。# 读取第一个波段并将其转为array band_1 = data.GetRasterBand(1) band_1 = band_1.ReadAsArray(0, 0, img_width, img_height) # 行，列2.2 写入影像

写一个遥感影像分为两部分：①完成这张图像，包括图像的行列号，像素点的数据类型。②影像在地球上的位置，即空间信息。包括仿射变换矩阵和地理位置。

# ①存储参数设置driver = gdal.GetDriverByName("GTiff")etype = gdal.GDT_Float32# ②orginal_file为原始影像，获取原始影像的空间参考和投影信息。将其写入新的影像proj = gdal.Open(orginal_file).GetProjection()transform = gdal.Open(orginal_file).GetGeoTransform()# 创建存储影像基本信息：# result_file结果影像# arr为需要存储的影像矩阵(行*列)，写如的格式为先列后行。band表示波段数量，etype是数据类型ds = driver.Create(result_file, arr.shape[1], arr.shape[0], band, etype) # eg.写单波段图像，并且将空间参考和投影信息写入。ds.GetRasterBand(1).WriteArray(arr)ds.SetGeoTransform(transform)ds.SetProjection(proj)ds.FlushCache()del ds3. 完整案例

两个函数，分别为文件的读写。此处随意选择一张tif格式的数据尝试，对需要更改的参数进行了详细的说明。

# 导入所需要的包from osgeo import gdalimport numpy as np# 读取遥感影像，获取影像的基本信息def read_tif(filename, b=1):'''filename为文件名称，b为波段数量；读取遥感影像并将其每个band转化为一列变量。b=1是默认参数，如果读取单波段影像，不需要输入；当输入影像为多波段时，需要将b改为影像的波段数。'''# 打开图像并创建空间 data = gdal.Open(filename) driver = gdal.GetDriverByName('GTiff') # driver = data.GetDriver() # 数据集的基本信息 print('Raster Driver : {d}\n'.format(d=driver.ShortName)) print('影像的波段数: ', data.RasterCount) img_width, img_height = data.RasterXSize, data.RasterYSize print('影像的列，行数: {r}rows * {c}colums'.format(r=img_width, c=img_height)) print('栅格数据的空间参考：{}'.format(data.GetGeoTransform())) # 栅格数据的6参数 print('投影信息：{}\n'.format(data.GetProjection())) # 栅格数据的投影 # 读取影像的元数据,获取各个波段的信息 print(data.GetMetadata()) band = None for i in range(b): band_1 = data.GetRasterBand(i + 1) band_1 = band_1.ReadAsArray(0, 0, img_width, img_height).flatten() # 行，列 if band is None: band = band_1 else: band = np.vstack((band, band_1)) return band.T def arr2raster(self, orginal_file, result_file, arr, band=1): ''' orginal_file为原始影像；目的为获取原始影像的地理参考。 result_file为要保存的文件的名称，arr为转存的矩阵，proj和transform分别为投影信息和仿射变换矩阵 band=1,表示默认情况生成单波段图像，若要多波段图像，在输入的时候更改band为要创建的波段数 ''' driver = gdal.GetDriverByName("GTiff") etype = gdal.GDT_Float32 proj = gdal.Open(orginal_file).GetProjection() transform = gdal.Open(orginal_file).GetGeoTransform() # no_data_value = '' if band == 1: ds = driver.Create(result_file, arr.shape[1], arr.shape[0], band, etype) # 行，列 # 写单波段图像 ds.GetRasterBand(1).WriteArray(arr) else: ds = driver.Create(result_file, arr.shape[2], arr.shape[1], band, etype) # 写多波段图像 for i in range(band): ds.GetRasterBand(i + 1).WriteArray(arr[i, :, :]) ds.SetGeoTransform(transform) ds.SetProjection(proj) ds.FlushCache() del ds# 案例演示：读取影像# 使用的是landsat5/TM影像，查看演示效果filename=r'C:\Users\lenovo\Desktop\test.tif'datasets = read_tif(filename, b=6)# 案例演示：写入影像# 实际中，会对影像进行操作，将结果再次转为tif。例如：将分类结果转为tif影像# arr为需要转存的矩阵，格式是行*列。此处将原始影像的band1转为tif，行列号721 * 1143arr = datasets[:,0].reshape(721, 1143)result_file = r'C:\Users\lenovo\Desktop\result.tif'arr2raster(filename, result_file, arr, band=1)