GDAL教程1
翻译+自己的理解
from __future__ import print_function
from osgeo import gdal
print("GDAL's version is:" + gdal.__version__)
print(gdal)
GDAL's version is:2.3.3
<module 'osgeo.gdal' from 'C:\\Users\\Dash\\Anaconda3\\lib\\site-packages\\osgeo\\gdal.py'>
我们使用上面这种办法来引入gdal模块
可以看到，我的GDAL版本是2.3.3
我们引入gdal子模块之后，我们就可以使用GDAL提供的API了
每次在GDAL中打开一个影像，我们就创建了一个GDALDataset对象。使用GDAL的Open函数去打开一张影像。
#Open a GDAL dataset
dataset = gdal.Open('tiffbigcut.tif', gdal.GA_ReadOnly)
print(dataset)
<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x0000019B43318E40> >
可以看到，已经有了一个dataset对象
下面查看一些图像的属性
num_bands = dataset.RasterCount#查看波段数
print("波段数为：{}\n".format(num_bands))#查看行列数
rows = dataset.RasterYSize
cols = dataset.RasterXSize
print('Image size is: {r} rows x {c} columns\n'.format(r=rows, c=cols))#查看描述和元数据
desc = dataset.GetDescription()#GetDescription方法
metadata = dataset.GetMetadata()#GetMetadata方法
print('Raster description: {desc}'.format(desc=desc))
print('Raster metadata:')
print(metadata)
print('\n')#查看打开这个影像的驱动
driver = dataset.GetDriver()#GetDriver方法
print('Raster driver: {d}\n'.format(d=driver.ShortName))#查看投影信息
proj = dataset.GetProjection()
print('Image projection:')
print(proj + "\n")#查看geo-transform
gt = dataset.GetGeoTransform()
print('Image geo-transform:{gt}\n'.format(gt=gt))
波段数为：3Image size is: 5140 rows x 5221 columnsRaster description: tiffbigcut.tif
Raster metadata:
{'AREA_OR_POINT': 'Area', 'TIFFTAG_XRESOLUTION': '1', 'TIFFTAG_YRESOLUTION': '1'}Raster driver: GTiffImage projection:
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],UNIT["degree",0.0174532925199433],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"]],AUTHORITY["EPSG","32650"]]Image geo-transform:(415875.0, 30.0, 0.0, 4064445.0, 0.0, -30.0)
上述属性的获取都是容易的。大家在使用的时候要注意大小写的问题哦（别问我怎么知道的）。
最后一个属性Geo-Transform地理变换。这组六个数字提供了转换像素和投影坐标所需要的信息。
GDALDataset对象包含了许多有用的信息，但是我们并不直接用他们来读取影像。我们需要去单独的去访问每个波段，利用的是GetRasterBand方法
blue = dataset.GetRasterBand(1)
print(blue)
<osgeo.gdal.Band; proxy of <Swig Object of type 'GDALRasterBandShadow *' at 0x0000019B45A03BA0> >
原图像是RGB的，1代表了blue。好，现在获得了blue，下面我们来看看blue的情况。
#查看该波段的数据类型
#返回一个数据类型的表示号码
datatype = blue.DataType
print('波段影像的类型是{dt}'.format(dt=datatype))#查看数据类型的名称
datatype_name = gdal.GetDataTypeName(blue.DataType)
print('波段类型名称是{}'.format(datatype_name))#产看该数据结构占用的存储空间
bytes = gdal.GetDataTypeSize(blue.DataType)
print('波段的数据结构的大小:{}'.format(bytes))#查看一些该波段的统计量
band_max, band_min, band_mean, band_stddev = blue.GetStatistics(0,1)
print('Band mean, stddev{m},{s}:\n'.format(m = band_mean, s=band_stddev))
波段影像的类型是2
波段类型名称是UInt16
波段的数据结构的大小:16
Band mean, stddev10242.75987232,1089.34635248:
对于大多数应用场景，我们需要使用gdal把我们的栅格数据读进内存。当我们装载我们的波段数据的时候，我们会将他读入numpy的二维数组。numpy的效率是很高的，能够提升我们的处理效率。
下面我们导入numpy
import numpy as np
print(np.__version__)
1.16.5
看到我的版本为1.16.5
使用ReadAsArray方法将数据读入numpy数组
blue_data = blue.ReadAsArray()
print(blue_data)
#查看一些属性
print('蓝波段的均值是:{m}'.format(m=blue_data.mean()))
print('尺寸是: {sz}'.format(sz = blue_data.shape))
[[ 9852 8887 9635 ... 9995 10235 10351][ 9894 9090 9692 ... 9978 10064 10066][ 9305 9039 9054 ... 9935 10079 10084]...[10260 10265 10308 ... 10097 11098 11570][10249 10237 10237 ... 10576 11155 11194][10224 10213 10187 ... 12177 11780 10591]]
蓝波段的均值是:10242.759872320477
尺寸是: (5140, 5221)
读入了一个numpy数组的数据，我们就可以通过Numpy来进行一些统计、以及线性代数、数组运算的工作。
下面，我们尝试将数据读入一个三维数组
#首先初始化一个三维的numpy数组
image = np.zeros((dataset.RasterYSize, dataset.RasterXSize, dataset.RasterCount))#好了，现在我们利用一个循环将数据读入我们初始化的numpy数组
for b in range(dataset.RasterCount):#注意，GDAL波段的索引从1开始！！！而Python是从0开始的！！！#注意上一行！！！！！！！！！！！！！！！band = dataset.GetRasterBand(b+1)#利用循环中的b来控制将波段数据读入numpy数组image[:, :, b] = band.ReadAsArray()#读完了，打印出来看看情况
#print(image)
#数据就不打印了，太长
print(image.dtype)
float64
我们现在做一个小小的调整，确保我们的GDAL影像能够匹配读入的Numpy数组的数据类型。GDAL恰有这样一个函数可以做GDAL到Numpy的转换，我们可以直接使用它而不用担心那些细节。
from osgeo import gdal_array#首先读取独立栅格波段数据的数据类型
image_datatype = dataset.GetRasterBand(1).DataType#下面进行numpy数组的初始化，但是这次初始化的时候，我们利用gdal的函数
#来确保数据类型一致
#这个函数就是gdal_array.GDALTypeCodeToNumericTypeCode（）
image_correct = np.zeros((dataset.RasterYSize,dataset.RasterXSize,dataset.RasterCount),dtype=gdal_array.GDALTypeCodeToNumericTypeCode(image_datatype))#好，此时再将数据读入数组，和上一次一样的
for b in range(dataset.RasterCount):band = dataset.GetRasterBand(b+1)image_correct[:, :, b] = band.ReadAsArray()print("对比数据类型：")
print("没有利用我们的这个函数的时候:{dt}".format(dt=image.dtype))
print("利用之后:{dt}".format(dt=image_correct.dtype))
对比数据类型：
没有利用我们的这个函数的时候:float64
利用之后:uint16
发现我们的numpy数组和原本的数据使用了相同的数据类型，因而通过这种方法，减少了4倍的内存空间！
我们下次将讨论如何取消内存分配，现在我们将他们设置为None
dataset = None
image = None
image_correct = None
翻译完了啊啊啊啊啊啊啊！