gdal库读取栅格影像数据

## 先介绍一下一些函数

gdal.Open('文件路径')       #切记gdal库的方法都是大写字母开头的

我们拿到dataset句柄之后，就可以对某个遥感影像进行具体操作了。

比如获取metadata信息，这个metadata信息其实我也不知道是怎么描述的，等到我之后学习了再来补吧。

dataset5.GetMetadata()  获取栅格数据元数据

from osgeo import gdal
import numpy as np
import cv2 as cv
dataset5=gdal.Open(r"D:\RS&GIS\landsat\LC81230372018098LGN00\LC08_L1TP_123037_20180408_20180417_01_T1_B5.TIF")
dataset4=gdal.Open(r'D:\RS&GIS\landsat\LC81230372018098LGN00\LC08_L1TP_123037_20180408_20180417_01_T1_B4.TIF')

print(dataset5.GetMetadata())

输出

{'AREA_OR_POINT': 'Point'}

• dataset5.RasterCount 获得栅格数据集的波段数

  print(dataset5.RasterCount)
  
  输出
  1
  

• dataset.GetGeoTransform() 栅格数据的六参数，其实就是地理仿射变换参数。在GDAL中，这六个参数包括左上角坐标，像元X、Y方向大小，旋转等信息。 要注意，Y方向的像元大小为负值。

• print(dataset5.GetGeoTransform())
  
  输出
  (158385.0, 30.0, 0.0, 3792015.0, 0.0, -30.0)

• GetProjection() 栅格数据的投影

• print(dataset5.GetProjection().strip().split(','))
  
  输出
  
  ['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"]]']
  
  进程已结束,退出代码0
  

gdalGetRasterBand(index)

这是获取遥感影像的某个波段的函数，和其他C数组类型不同，里面的index是从1开始的。这里来演示一下，我这里是直接打开两个单波段的影像。

from osgeo import gdal
import numpy as np
import cv2 as cv
dataset5=gdal.Open(r"D:\RS&GIS\landsat\LC81230372018098LGN00\LC08_L1TP_123037_20180408_20180417_01_T1_B5.TIF")

data5=dataset5.GetRasterBand(1)

输出
D:\desktop\pythonProject1\venv\Scripts\python.exe D:/desktop/pythonProject1/程序/地理数据简单操作.py

进程已结束,退出代码0

这里我们获取某个波段之后，我们还可以做什么呢？

当然是得到某个波段的像元信息呀

GDAL提供了下面两个函数来访问影像的数值。

ReadRaster() 读取图像数据(以二进制的形式)

ReadAsArray() 读取图像数据(以数组的形式)根据个人经验，一般是使用ReadAsArray(),这可以和numpy库通过数值来进行连接，这个函数很神奇可以当做一个伸缩自如的窗口，但是之前你得给它一个初始位置

ReadRaster(xoff=0, yoff=0, xsize=None, ysize=None, buf_xsize=None, buf_ysize=None, buf_type=None,
 band_list=None, buf_pixel_space=None, buf_line_space=None, 
buf_band_space=None, resample_alg=0, 
callback=None, callback_data=None) method of osgeo.gdal.Dataset instance

参数具体解释我从官网上直接摘抄下来了。

• xoff,yoff ：指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位）。

• xsize,ysize ： 指定要读取部分图像的矩形的长和宽（以像元为单位）。

• buf_xsize，buf_ysize ：可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高，gdal将帮你缩放到这个大小。

• buf_type ：可以对读出的数据的类型进行转换（比如原图数据类型是short，你要把它们缩小成byte）。

• band_list ：适应多波段的情况。可以指定要读取的波段。

下面来简单演示一下，这是一幅原始没有结果任何处理的影像，专业人士一眼就可以看出来了

​
dn5=data5.ReadAsArray(2000,2000,5,5,10,10)
print(dn5)

输出
[[21244 21244 21172 21172 23818 23818 24048 24048 24507 24507]
 [21244 21244 21172 21172 23818 23818 24048 24048 24507 24507]
 [16414 16414 16898 16898 22044 22044 24571 24571 24690 24690]
 [16414 16414 16898 16898 22044 22044 24571 24571 24690 24690]
 [16079 16079 16758 16758 21536 21536 24756 24756 24973 24973]
 [16079 16079 16758 16758 21536 21536 24756 24756 24973 24973]
 [16497 16497 16863 16863 22777 22777 25708 25708 25722 25722]
 [16497 16497 16863 16863 22777 22777 25708 25708 25722 25722]
 [18269 18269 18245 18245 23081 23081 24775 24775 24815 24815]
 [18269 18269 18245 18245 23081 23081 24775 24775 24815 24815]]

​

大家想一下，万一我不是从第1200行，1200列开始呢，要是从一个很大的行列号开始的话肯定会超出数据范围。