1. gdal介绍

GDAL(Geospatial Data Abstraction Library)主要用来读取地理空间数据    
    
    
，现在的GDAL包并不是单独的GDAL 


     


     


，而是集成了GDAL和OGR的
    
    
。OGR用于处理矢量数据     


     


     。因此

 




 




 ，GDAL既可以用来处理栅格也可以处理矢量文件 




 




 



。

2. 代码详解

下面代码主要为读写tif的代码   

   

   
，图像其实就是一个二维矩阵  

     

     

，其他格式的遥感影像都可以参考

   

   

。包括nc


  




  




 ，hdf等等     

     

    。

2.1 读取数据

# 调包 from osgeo import gdal import 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 ds

3. 完整案例

两个函数










，分别为文件的读写

   

   

。此处随意选择一张tif格式的数据尝试               ，对需要更改的参数进行了详细的说明     

     

    。

# 导入所需要的包 from osgeo import gdal import 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=rC:\Users\lenovo\Desktop\test.tif datasets = read_tif(filename, b=6) # 案例演示：写入影像 # 实际中




    




    




 ，会对影像进行操作










，将结果再次转为tif   




   




   

。例如：将分类结果转为tif影像 # arr为需要转存的矩阵               ，格式是行*列



 



 



。此处将原始影像的band1转为tif



     



     



 ，行列号721 * 1143 arr = datasets[:,0].reshape(721, 1143) result_file = rC:\Users\lenovo\Desktop\result.tif arr2raster(filename, result_file, arr, band=1)