Python——使用GDAL操作栅格数据（自用）

 range()函数
range()函数返回的是可迭代对象，不是列表类型

val1=range(5,0,-1)
print(val1)

(1) range(stop)，例如range(10)是默认从0开始，步长为1，不包括10

(2) range(start,stop) ,其运行结果则由起始值开始，步长为1，不包括停止值

(3)range(start,stop,step),创建一个在[start,stop)之间，步长为step

tostring()方法可以将NumPy数组转化为二进制p66

读取遥感影像信息

>>> from osgeo import gdal
>>> dataset = gdal.Open("F:/data/BIL_BANKUAI.tif")
>>> dir(dataset)  #该tif文件可执行的操作
['AbortSQL', 'AddBand', 'AddFieldDomain', 'AdviseRead', 'BeginAsyncReader', 'BuildOverviews', 'ClearStatistics', 'CommitTransaction', 'CopyLayer', 'CreateLayer', 'CreateMaskBand', 'DeleteLayer', 'EndAsyncReader', 'ExecuteSQL', 'FlushCache', 'GetDescription', 'GetDriver', 'GetFieldDomain', 'GetFileList', 'GetGCPCount', 'GetGCPProjection', 'GetGCPSpatialRef', 'GetGCPs', 'GetGeoTransform', 'GetLayer', 'GetLayerByIndex', 'GetLayerByName', 'GetLayerCount', 'GetMetadata', 'GetMetadataDomainList', 'GetMetadataItem', 'GetMetadata_Dict', 'GetMetadata_List', 'GetNextFeature', 'GetProjection', 'GetProjectionRef', 'GetRasterBand', 'GetRootGroup', 'GetSpatialRef', 'GetStyleTable', 'GetSubDatasets', 'GetTiledVirtualMem', 'GetTiledVirtualMemArray', 'GetVirtualMem', 'GetVirtualMemArray', 'IsLayerPrivate', 'RasterCount', 'RasterXSize', 'RasterYSize', 'ReadAsArray', 'ReadRaster', 'ReadRaster1', 'ReleaseResultSet', 'ResetReading', 'RollbackTransaction', 'SetDescription', 'SetGCPs', 'SetGeoTransform', 'SetMetadata', 'SetMetadataItem', 'SetProjection', 'SetSpatialRef', 'SetStyleTable', 'StartTransaction', 'TestCapability', 'WriteArray', 'WriteRaster', '_SetGCPs', '_SetGCPs2', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__swig_destroy__', '__weakref__', 'this', 'thisown']
>>> dataset.GetRasterBand(1) #查看波段
<osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x0000018FB73A1360> >
>>> dataset.GetMetadata() #访问元数据
{'AREA_OR_POINT': 'Area', 'TIFFTAG_XRESOLUTION': '1', 'TIFFTAG_YRESOLUTION': '1'}
>>> dataset.GetDescription() #获取栅格数据信息
'F:/data/BIL_BANKUAI.tif'
>>> dataset.RasterCount #波段数查询
1
>>> dataset.RasterXSize,dataset.RasterYSize  #获取影像行列数目
(1784, 1844)
>>> band.DataType #
1

创建与保存栅格数据集

1.使用CreateCopy()方法创建影像

2. 使用Create()方法创建影像

>>> from osgeo import gdal
#利用Create()函数创建GeoTIFF格式单波段栅格影像，大小512*512像元，数据类型为Byte
>>> driver=gdal.GetDriverByName('GTiff')
>>> dst_filename='E:/data/x_tmp.tif'
>>> dst_ds=driver.Create(dst_filename,512,512,1,gdal.GDT_Byte)

#对创建的单波段影像的空间范围与坐标系进行设置
>>> import numpy
>>> from osgeo import osr
>>> dst_ds.SetGeoTransform([444720,30,0,3751320,0,-30])
0
>>> srs=osr.SpatialReference()
>>> srs.SetUTM(11,1)
0
>>> srs.SetWellKnownGeogCS('NAD27')
0
>>> dst_ds.SetProjection(srs.ExportToWkt())
0
>>> raster=numpy.zeros((512,512))
>>> dst_ds.GetRasterBand(1).WriteArray(raster)

3.创建多波段影像

#建立3波段GeoTIFF
>>> from osgeo import gdal
>>> import numpy
>>> dataset=gdal.Open("E:/data/BIL_BANKUAI.tif")
>>> width=dataset.RasterXSize
>>> height=dataset.RasterYSize  

#出现错误ImportError: numpy.core.multiarray failed to import，则需要将numpy卸载重新安装，或升级其版本
>>> datas=dataset.ReadAsArray(0,0,width,height)
>>> driver=gdal.GetDriverByName("GTiff")
>>> tods=driver.Create("E:/data/geotiff_BIL_BANKUAI_3.tif",width,height,3,options=["INTERLEAVE=PIXEL"])
#band_list参数是由波段数决定的，需根据源数据进行判断
#>>> tods.WriteRaster(0,0,width,height,datas.tostring(),width,height,band_list=[3,4,5])
#ERROR 1: Buffer too small
#3
>>> tods.WriteRaster(0,0,width,height,datas.tostring(),width,height,band_list=[3])
0
>>>tods.FlushCache()

4.GDAL写操作时的空间投影处理

#使用GDAL创建数据时影像的空间参数信息是单独处理的
#导入NAD27空间参考,NAD27是美国查工空间参考系统代码
>>> from osgeo import osr
>>> from osgeo import gdal
>>> dataset=gdal.Open("E:/data/BIL_BANKUAI.tif")
>>> width=dataset.RasterXSize
>>> height=dataset.RasterYSize
>>> datas=dataset.ReadAsArray(0,0,width,height)
>>> driver=gdal.GetDriverByName("GTiff")
>>> tods=driver.Create("E:/data/x_BIL_BANKUAI_3.tif",width,height,3,options=["INTERLEAVE=PIXEL"])
#（左上角x坐标， 水平分辨率，旋转参数， 左上角y坐标，旋转参数，-垂直分辨率）
>>> tods.SetGeoTransform([444720,30,0,3751320,0,-30])
0
>>> srs=osr.SpatialReference()
>>> srs.SetUTM(11,1)
0
>>> srs.SetWellKnownGeogCS('NAD27')
0
>>> tods.SetProjection(srs.ExportToWkt())
0
#若使用空间参考比较麻烦，则可只定义六参数变换

5.建立影像金字塔

#建立ERDAS Images风格影像金字塔

#设置影像金字塔风格
>>> from osgeo import gdal
>>> gdal.SetConfigOption('HFA_USE_RRD','YES')

#建立影像金字塔
>>> dataset=gdal.Open("E:/data/BIL_BANKUAI.tif")
#list里面是要生成金字塔的层级，全部为2的幂次，第一层是原始数据
>>> dataset.BuildOverviews(overviewlist=[2,4,8,16,32,64,128])
0

GDAL的其他操作

1.分别利用GDAL和Pillow读取数据

#利用Pillow读取影像数据
>>> from PIL import Image
>>> im=Image.open("E:/data/BIL_BANKUAI.tif")  #im可类比为GDAL的数据集
>>> region=im.crop((30,70,35,75))
>>> region.tobytes()
b'\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01'
>>>

2.Pillow与GDAL读取数据的转换

#GDAL与Pillow的转换
>>> import numpy as np
>>> data=dataset.ReadAsArray(30,70,5,5)   #GDAL的参数为（原点X，原点Y，宽，高）
>>> from numpy import reshape
>>> datas=[reshape(i,(-1,1)) for i in data]
>>> datas=np.concatenate(datas,1)
>>> datas.tostring()
b'\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01\x01'

3.ColorMap颜色定义

>>> from osgeo import gdal
>>> from osgeo import gdalconst
>>> dataset=gdal.Open('E:/data/BIL_BANKUAI.tif')
>>> colormap=dataset.GetRasterBand(1).GetRasterColorTable()   #GetRasterColorTable()获得颜色表
>>> dir(colormap)
['Clone', 'CreateColorRamp', 'GetColorEntry', 'GetColorEntryAsRGB', 'GetCount', 'GetPaletteInterpretation', 'SetColorEntry', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__swig_destroy__', '__weakref__', 'this', 'thisown']
>>> colormap.GetCount()  #GetCount获得颜色数量
256
>>> colormap.GetPaletteInterpretation()   #从GetPaletteInterpretation()可知该影像颜色表为一个RGB颜色表
1
>>> gdalconst.GPI_Gray  
0                       #调色板类型值，Gray灰度图像对应0，RGB（RBGA图像）对应1，CMYK图像对应2，HLS图像对应3
>>> gdalconst.GPI_RGB
1

#颜色查找表进行遍历操作
>>> for i in range(colormap.GetCount()):
...     print("{}:{}".format(i,colormap.GetColorEntry(i)))
...
0:(0, 0, 0, 255)      #四个值在不同的调色板中具有不同的定义，见p77表2.6
1:(0, 138, 0, 255)
2:(237, 153, 0, 255)
3:(0, 254, 0, 255)
4:(126, 254, 211, 255)
5:(254, 0, 0, 255)
6:(0, 0, 0, 255)
7:(0, 0, 0, 255)
8:(0, 0, 0, 255)
9:(0, 0, 0, 255)
10:(0, 0, 0, 255)
11:(0, 0, 0, 255)
...

4.地图代数的计算

#以计算NDVI为例

#将波段2读入数组data2，将波段3读入数组data3
>>> from osgeo import gdal
>>> import numpy
>>> dataset=gdal.Open('E:/data/BIL_BANKUAI.tif')
#>>> band=dataset.GetRasterBand(1)
>>> band2=dataset.GetRasterBand(2)
>>> band3=dataset.GetRasterBand(3)
>>> cols=100
>>> rows=100
#ReadAsArray(原点X,原点Y，宽，高)
>>> data2=band2.ReadAsArray(0,0,cols,rows).astype(numpy.float16)
>>> data3=band3.ReadAsArray(0,0,cols,rows).astype(numpy.float16)
#去除data2和data3中的0值，以防出现错误导致崩溃,Greater函数用来判断参数一是否大于参数二
>>> mask=numpy.greater(data2+data3,0)
#计算NDVI
#choose(a, choices, out=None, mode='raise')
#参数 a ：它必须是一个 int 型的 数组，并且 a 中的元素，必须是0~n-1之间的数，这里的n表示的就是组choices数组最外层的维度数。
#choices：表示的是要操作的数组，要注意的是choices的数组的维度是一定要和a进行匹配的，如果匹配不了，会出现错误。
#参数out：接收运算结果的数组，它的维度一定要和 a 是一样的，是可选参数。
#参数mode：默认的是raise，表示的是a数组中的元素不能超过 n ，她还有两个可选值；clip：将 a 中的 元素 如果小于0，则将其变为0，如果大于n-1，则变为n-1；wrap：将a中的值 value变为value mod n，即值除以n的余数。
>>> ndvi=numpy.choose(mask,(-99,(data3-data2)/(data3+data2)))
>>> ndvi
array([...])