gdal概览

1 gdal库

2 栅格驱动

读取栅格数据的流程：
（1）获取栅格驱动(会有1个默认的，可以不创建，最好创建上，看着清晰)；
（2）驱动进行注册（会有默认值，可以不创建，创建的话可以选择只读，或者读写都可以，相当于给栅格驱动设置一个权限）；
（3）通过驱动获取数据集；
（4）通过数据集操作栅格属性；

try:
    from osgeo import gdal,ogr   # gdal是处理栅格，ogr处理矢量
except:
    import gdal

driverCount = gdal.GetDriverCount()
print(driverCount)

#gdal.AllRegister()
driver = gdal.GetDriverByName('GTiff')
driver.Register()

print(driver.ShortName)
print(driver.LongName)

print(dir(driver))

3 栅格数据集（就是包含各种栅格属性的一个类）

3.1 坐标（6个参数）

# 读取地理坐标
from osgeo import gdal

filePath = '1.tif'  # tif文件路径
dataset = gdal.Open(filePath)  # 打开tif

adfGeoTransform = dataset.GetGeoTransform()  # 读取地理信息

# 左上角地理坐标
print(adfGeoTransform[0])
print(adfGeoTransform[3])

nXSize = dataset.RasterXSize  # 列数
nYSize = dataset.RasterYSize  # 行数

print(nXSize, nYSize)

arrSlope = []  # 用于存储每个像素的（X，Y）坐标
for i in range(nYSize):
    row = []
    for j in range(nXSize):
        px = adfGeoTransform[0] + i * adfGeoTransform[1] + j * adfGeoTransform[2]
        py = adfGeoTransform[3] + i * adfGeoTransform[4] + j * adfGeoTransform[5]
        col = [px, py]  # 每个像素的经纬度
        row.append(col)
        print(col)
    arrSlope.append(row)

上面的代码其实已经实现获取tif中经纬度，如果大家仔细研究一下会发现，其实我们使用的就是gdal里面的GetGeoTransform方法读取坐标，简单介绍一下该方法，该方法会返回以下六个参数

GT(0) 左上像素左上角的x坐标。
GT(1) w-e像素分辨率/像素宽度。
GT(2) 行旋转（通常为零）。
GT(3) 左上像素左上角的y坐标。
GT(4) 列旋转（通常为零）。
GT(5) n-s像素分辨率/像素高度（北上图像为负值）

3.1.2 tif文件的地理坐标（两种情况）

一是只有tif文件，地理坐标存储在这个tif文件里；
二是有附带的tfw文件，这个文件里存储的地理坐标的6个参数；
arcgis读取tif文件的时候首先看tif文件里有没有地理坐标，有的话直接用，没有的话才会去找tfw文件；

3.2 波段数、大小、投影等信息

try:
    from osgeo import gdal,ogr
except:
    import gdal

gdal.AllRegister()
dataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif')

print(dir(dataset))

#波段信息
bandCount = dataset.RasterCount
print(bandCount)

#大小
weight,height = dataset.RasterXSize,dataset.RasterYSize
print(weight,height)

#空间信息
transform = dataset.GetGeoTransform()
print(transform)

#投影信息
proj = dataset.GetProjection()

#描述信息
descrip = dataset.GetDescription()
print(descrip)

3.3 读取栅格像元

读取栅格数据的流程：
（1）获取栅格驱动(会有1个默认的，可以不创建)；
（2）驱动进行注册（会有默认值，可以不创建，创建的话可以选择只读，或者读写都可以，相当于给栅格驱动设置一个权限）；
（3）通过驱动获取数据集；
（4）通过数据集操作栅格属性；

try:
    from osgeo import gdal,ogr
except:
    import gdal

gdal.AllRegister()      # 只读，这里前面没有指定驱动，它自动会用tif
dataset = gdal.Open('C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif')

width = dataset.RasterXSize       # 数据集中的各种属性信息
height = dataset.RasterYSize
bandCont = dataset.RasterCount

for i in range(bandCont):
    band = dataset.GetRasterBand(i+1)      # 获取波段，波段从1开始，不是从0开始
    bandinform = band.ReadRaster(0,0,3,3)  # 获取栅格数值0-3
    print(bandinform)
    #获得波段的数值类型
    dataType = band.DataType
    print(dataType)

    #获得影像的nodata
    nodata = band.GetNoDataValue()
    print(nodata)

    #获得最大值最小值，这个波段中的最大最小值
    maxmin = band.ComputeRasterMinMax()
    print(maxmin)

    imageArray = band.ReadAsArray(0,0,3,3,10,10)   # 读取栅格像元并输出为numpy的array形式
    print(imageArray)

3.4 创建栅格影像

创建栅格影像流程：
（1）创建驱动；
（2）定义数据集名字，数据集影像宽和高；
（3）创建数据集；
（4）添加栅格像元值；


上面两种读写方式在下面的例子中都有实际运用：
（1）WriteRaster: 看3.4.4 随机裁剪栅格；
（2）WriteArray: 看3.4.5 计算NDVI波段。

3.4.1 直接用数组创建数据集

像元值存储的类型，不同数值类型得到的像元值范围不一样。

from osgeo import gdal
import numpy as np

srcFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1.tif'
dstFile = 'C:\\Users\\Administrator\\Desktop\\gdal\\样例数据\\GF1\\GF1Create.tif'
dataSet = gdal.Open(srcFile)

width = dataSet.RasterXSize
height = dataSet.RasterYSize
bandCount = dataSet.RasterCount

driver = gdal.GetDriverByName('GTiff')
print(dir(driver))
metadata = driver.GetMetadata()
if gdal.DCAP_CREATE in metadata and metadata['DCAP_CREATE'] == 'YES':
    print('支持create方法')
else:
    print('不支持create方法')
if gdal.DCAP_CREATECOPY in metadata and metadata['DCAP_CREATECOPY'] == 'YES':
    print('支持createCopy方法')
else:
    print('不支持createCopy方法')
dataSetOut = driver.Create(dstFile,width,height,bandCount,gdal.GDT_Byte)   # 8位无符号整数

#空间参考
srcProj