GDAL的使用

• 栅格位置(像素或者是行坐标)和地理参考坐标之间的转换可以通过仿射变换实现，仿射矩阵可以通过GDALDataset::GetGeoTransform()得到，依据下面的公式将像素/行坐标转换到地理参考空间：
  $$\begin{gathered} X_{geo}=GT(0)+Xpixel.GT(1)+Yline.GT(2) \\ {Y}_{geo}=GT(3)+Xpixel.GT(4)+Yline.GT(5) \end{gathered}$$
  其中像素坐标左上角为$(0,0)$，$(GT(0),GT(3))$是栅格左上角坐标，$GT(1)$是像素宽度，$GT(5)$是像素的高度；
• 可以通过地理控制点(GCP)来描述栅格数据集地理参考，关联了栅格位置和地理参考系统的一个或者多个位置，一个数据集会有一套控制点，控制点通过GDALDataset::GetGCPProjection()得到，每个控制点包含Char pszId(控制点标识符)，Char pszInfo(通常是空串)，double dfGCPPixel+double dfGCPLine(像素、线的位置是控制点在栅格的位置)，double dfGCPX+double dfGCPY+double dfGCPZ(地理参考位置)，本条和上一条都描述了栅格位置和地理参考坐标之间的关系；
• 栅格波段GDALRasterBand
• 颜色表：像素值被用作颜色表的下标；

try:
    import gdal
except:
    from osgeo import gdal	# gdal1.6
from osgeo.gdalconst import *	# 常用的常量

# 要读取数据之前需要载入数据驱动
# 注册所有数据驱动
gdal.AllRegister()
# 某一类型的数据驱动，按照数据格式加载
driver = gdal.GetDriverByName('GTiff')
driver.Register()

# 查看系统支持的数据格式，其中shortname就是上面GetDriverByName的参数，对于不同的gdal版本GetDriver的结果可能不同
drv_count = gdal.GetDriverCount()
for idx in range(drv_count):
		driver = gdal.GetDriver(idx)
		print( "%10s: %s" % (driver.ShortName, driver.LongName))

# 读取遥感影像
# 打开GeoTIF文件
dataset = gdal.Open("/gdata/geotiff_file.tif")
# python的dir()函数可以快速查看某对象可用的操作
dataset.GetDescription() 	# 获得栅格的描述信息
dataset.RasterCount 	# 获得栅格数据集的波段数
dataset.RasterXSize 	# 栅格数据的宽度(X方向上的像素个数)
dataset.RasterYSize 	# 栅格数据的高度(Y方向上的像素个数)
dataset.GetGeoTransform() 	# 栅格数据的六参数，这六个参数包括 左上角坐标 ， 像元X、Y方向大小 ， 旋转 等信息。 要注意， Y 方向的像元大小为负值
GetProjection() 	# 栅格数据的投影
dataset.GetMetadata()		# 获取元数据

# 获取数据集的信息
band = dataset.GetRasterBand(1)		# 获取第一个波段，下标从0开始
band.XSize, band.YSize, band.DataType	# 宽高和数据类型
# 数据类型对应
# 未知或未指定类型 gdalconst.GDT_Unknown 0
# 8位无符整型 gdalconst.GDT_Byte 1
# 16位无符整型 gdalconst.GDT_UInt16 2
# 16位整型 gdalconst.GDT_Int16 3
# 32位无符整型 gdalconst.GDT_UInt32 4
# 32位整型值 gdalconst.GDT_Int32 5
# 32位浮点型 gdalconst.GDT_Float32 6
# 64位浮点型 gdalconst.GDT_Float64 7
# 16位复数整型 gdalconst.GDT_CInt16 8
# 32位复数整型 gdalconst.GDT_CInt32 9
# 32位复数浮点型 gdalconst.GDT_CFloat32 10
# 64位复数浮点型 gdalconst.GDT_CFloat64 11
band.GetNoDataValue()		# 无意义填充值
band.GetMaximum()
band.GetMinimum()
band.ComputeRasterMinMax()

# 访问栅格数据集的数据
dataset.ReadRaster() 	# 读取图像数据(以二进制的形式)
dataset.ReadAsArray() 	# 读取图像数据(以数组的形式)，返回的是numpy的array
# 参数
# xoff,yoff ：指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位）
# xsize,ysize ： 指定要读取部分图像的矩形的长和宽（以像元为单位）
# buf_xsize，buf_ysize ：可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高， GDAL 将帮你缩放到这个大小
# buf_type ：可以对读出的数据的类型进行转换（比如原图数据类型是short，你要把它们缩小成byte）
# band_list ：适应多波段的情况。可以指定要读取的波段

# 查看图片信息
!gdalinfo /gdata/lu75i1.tif	
# 访问索引图像：所读的数据知识真实数据的索引，而不是灰度图像
dataset = gdal.Open('/gdata/lu75i1.tif')
band = dataset.GetRasterBand(1)
band.GetRasterColorInterpretation()		# 返回2，gdalconst.GCI_PaletteIndex，表示索引图
colormap = band.GetRasterColorTable()	# 获取颜色表
colormap.GetPaletteInterpretation()		# 获取颜色表的类型
colormap.GetCount()	# 颜色数量
for i in range(colormap.GetCount() - 10, colormap.GetCount()):
   print("%i:%s" % (i, colormap.GetColorEntry(i)))	# 获得颜色的四值元祖，例如rgb,cmyk
# 我们通过ReadRaster读出的数据值只是对应到这个表的一个索引而已。 
# 我们需要通过读出这些数据，并在真实数据表中找出真实数据， 重新组织成一个RGB表才能用来绘制。
# 如果不经过对应， 绘制出来的东西可能没有任何意义
# GTiff颜色表存储时是16位的，但是读取之后自动进行了处理变为0-255

# 创建影像
driver = gdal.GetDriverByName( 'GTiff' )
dst_filename = '/tmp/x_tmp.tif'
dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )
from osgeo import osr
dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] )
srs = osr.SpatialReference()
srs.SetUTM( 11, 1 )
srs.SetWellKnownGeogCS( 'NAD27' )
dst_ds.SetProjection( srs.ExportToWkt() )
raster = numpy.zeros( (512, 512) )
dst_ds.GetRasterBand(1).WriteArray( raster )

ref:
https://www.osgeo.cn/pygis/gdal.html