我要投稿
投稿指南
RSS订阅
3. 快速开始
其实在主站的教程里已经有python的示例了。但是我们还是按照自己的思路来开始吧。
第一步就是打开一个数据集。对于“数据集”这个名词大家可能不会太习惯，但是对于一般的格式来说，一个“数据集”就是一个文件，比如一个gif文件就是一个以gif为扩展名的文件。但是对于众多RS数据来说，一个数据集包含的绝对不仅仅是一个文件。对于很多RS数据，他们把一张图像分成数个图像文件，然后放在一个文件夹中，用一些额外的文件来组织它们之间的关系，形成一个“数据集”。如果你不理解，那么就算了，当成jpg或者gif文件好了。
下面我们打开一个tiff文件(GeoTIFF)。这个文件是我从GRASS的示例数据spearfish中导出的一个同名影像数据。
1 >>> import gdal
2 >>> dataset = gdal.Open("j:/gisdata/gtif/spot.tif")
3 >>> dir(dataset)
4 ['AddBand', 'AdviseRead', 'BuildOverviews', 'FlushCache', 'GetDescription', 'Get
5 Driver', 'GetGCPCount', 'GetGCPProjection', 'GetGCPs', 'GetGeoTransform', 'GetMe
6 tadata', 'GetProjection', 'GetProjectionRef', 'GetRasterBand', 'GetSubDatasets',
7 'RasterCount', 'RasterXSize', 'RasterYSize', 'ReadAsArray', 'ReadRaster', 'Refr
8 eshBandInfo', 'SetDescription', 'SetGCPs', 'SetGeoTransform', 'SetMetadata', 'Se
9 tProjection', 'WriteRaster', '__del__', '__doc__', '__init__', '__module__', '_b
10 and', '_o']
11 >>>
这样我们就打开了这个文件。并且我们可以看到可以供我们调用的函数们(更具体的API列表可以看这里)。现在我们不做修改，不做添加，所以只要带有Set开头的函数以及有Write开头的函数我们暂时都不管。因为RS影像必然要和地理上的位置挂上钩，才能把图像正确铺展到一个坐标系中。其中的信息和对应关系有点复杂，不适合在快速开始中介绍，我们暂时也先不管。这里需要注意的就是几个函数。
GetDescription 获得栅格的描述信息。
1 >>> dataset.GetDescription()
2 'j:/gisdata/gtif/spot.tif'
3 >>>
看来这里的图像描述是图像的路径名，但是这是和各种不同数据集相关的，不同数据集可能有不同的描述。这要看读取驱动的实现作者的高兴了。
RasterCount 获得栅格数据集的波段数。
GetRasterBand 获得栅格数据集的波段。
1 >>> dataset.RasterCount
2 1
3 >>> band = dataset.GetRasterBand(1)
4 >>>
这里需要解释的是Band这个词。这个词可以翻译成“波段”，“通道”等等。我这里把它统一称为“波段”。因为遥感卫星的传感器有很多个。一个传感器只负责接收一个频率范围的地物反射光波，一个频率范围的光波记录称为一个波段。是不是晕了？其实说得简单一点。其实你可以把波段看成红绿蓝几种颜色。图像不是分RGB三色吗？把R，G，B值都提取出来成为三个表。R值表就是波段一，G值表就是波段二，B值表就是波段三。
这里我们看到这张图只有一个波段(一种颜色)。就可以把它看成是一个灰度图(类似黑白照片)。如果RasterCount是3，就有可能是彩色图。如果RasterCount是比3大的数，恭喜你，你看到一张遥感影像。有很多卫星的传感器大于3个，比如TM就有7个波段，不仅有可见光，还有红外等其他非可见光。，所以，波段一般比RGB能表达的丰富地多。不过这样一来就需要我们从中挑出3个波段然后组合成RGB，当然这样就有可能使图像显示出来的东西不像平常我们看到的那样。这样安排是因为对科学有帮助(一些波段在科学家眼里比真实的彩色照片有价值)。不理解就跳过，很正常，我第一次听这种东西也觉得很玄:)
这里我们获取了第一个波段(红色值组成的表)。注意！这里的波段获取和通常的C数组获取不一样，开始是1不是0。获取了波段，我们就可以在下面的操作中读取这个波段的所有数值。
RasterXSize 图像的宽度(X方向上的像素个数)
RasterYSize 图像的高度(Y方向上的像素个数)
1 >>> dataset.RasterXSize
2 950
3 >>> dataset.RasterYSize
4 700
5 >>>
可以看出我们的图像大小是950*700。还是很小的一张图。
ReadRaster 读取图像数据(以二进制的形式)
ReadAsArray 读取图像数据(以数组的形式)
1 >>> help(dataset.ReadRaster)
2 Help on method ReadRaster in module gdal:
3 ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t
4 ype=None, band_list=None) method of gdal.Dataset instance
5 >>> help(dataset.ReadAsArray)
6 Help on method ReadAsArray in module gdal:
7 ReadAsArray(self, xoff=0, yoff=0, xsize=None, ysize=None) method of gdal.Dataset
8 instance
9 >>>
这两个函数很重要，它们直接读取图像的数据，可以看到两个函数的帮助中有一大溜的参数。解释一下： xoff,yoff,xsize,ysize 你可能不想读取整张图像。只想读取其中的一部分。那么就用xoff，yoff指定想要读取的部分原点位置在整张图像中距离全图原点的位置。用xsize和ysize指定要读取部分图像的矩形大小。
buf_xsize buf_ysize 你可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高，gdal将帮你缩放到这个大小。
buf_type 如果你要读取的图像的数据类型不是你想要的(比如原图数据类型是short，你要把它们缩小成byte)，就可以设置它。
band_list 这就适应上面多波段的情况。你可以指定读取的波段序列。要哪几个波段，不要哪几个波段，你说了算。
举个例子吧：
1 >>> dataset.ReadAsArray(230,270,10,10)
2 array([[255, 255, 255, 232, 232, 255, 255, 255, 255, 222],
3 [255, 255, 255, 255, 255, 255, 210, 110, 11, 122],
4 [255, 255, 255, 255, 255, 255, 210, 255, 11, 243],
5 [201, 255, 255, 255, 255, 200, 200, 110, 122, 243],
6 [111, 211, 255, 201, 255, 255, 100, 11, 132, 243],
7 [255, 100, 100, 100, 110, 100, 110, 111, 122, 243],
8 [255, 255, 255, 255, 255, 255, 122, 222, 255, 255],
9 [255, 255, 255, 255, 255, 255, 243, 243, 255, 255],
10 [255, 255, 255, 255, 255, 255, 255, 255, 255, 255],
11 [255, 255, 255, 255, 255, 255, 255, 255, 255, 255]],'b')
12 >>> dataset.ReadRaster(230,270,10,10)
13 '\xff\xff\xff\xe8\xe8\xff\xff\xff\xff\xde\xff\xff\xff\xff\xff\xff\xd2n\x0bz\xff\
14 xff\xff\xff\xff\xff\xd2\xff\x0b\xf3\xc9\xff\xff\xff\xff\xc8\xc8nz\xf3o\xd3\xff\x
15 c9\xff\xffd\x0b\x84\xf3\xffdddndnoz\xf3\xff\xff\xff\xff\xff\xffz\xde\xff\xff\xff
16 \xff\xff\xff\xff\xff\xf3\xf3\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff
17 \xff\xff\xff\xff\xff\xff\xff\xff\xff'
18 >>>
我们就把图像中位于230,270，宽度10高度10的数据读取出来了。
我们看完了数据集的主要函数。似乎已经够用了。的确，如果只是为了显示图像，这些的确已经够了。但是如果需要更多信息，我们就不得不进入波段操作数据(实际上我们大多数时候都需要进入band获取信息)。下面我们现在来看看刚才读取出来的那个band有些什么东西可以供我们操作的(具体的API列表看这里)。
1 >>> dir(band)
2 ['AdviseRead', 'Checksum', 'ComputeBandStats', 'ComputeRasterMinMax', 'DataType'
3 , 'Fill', 'FlushCache', 'GetDefaultHistogram', 'GetDescription', 'GetHistogram',
4 'GetMaximum', 'GetMetadata', 'GetMinimum', 'GetNoDataValue', 'GetOffset', 'GetO
5 verview', 'GetOverviewCount', 'GetRasterColorInterpretation', 'GetRasterColorTab
6 le', 'GetScale', 'GetStatistics', 'ReadAsArray', 'ReadRaster', 'SetDefaultHistog
7 ram', 'SetDescription', 'SetMetadata', 'SetNoDataValue', 'SetRasterColorInterpre
8 tation', 'SetRasterColorTable', 'WriteArray', 'WriteRaster', 'XSize', 'YSize', '
9 __doc__', '__init__', '__module__', '_o']
10 >>>
挑几个有用的吧。
1 >>> band.XSize
2 950
3 >>> band.YSize
4 700
5 >>> band.DataType
6 1
7 >>>
不用解释了吧，波段图像的宽和高(象元为单位)。DataType，图像中实际数值的数据类型。具体数据类型定义在gdalconst模块里。使用的时候用import gdalconst引入。
1 >>> import gdalconst
2 >>> dir(gdalconst)
3 ['CE_Debug', 'CE_Failure', 'CE_Fatal', 'CE_None', 'CE_Warning', 'CPLES_Backslash
4 Quotable', 'CPLES_CSV', 'CPLES_SQL', 'CPLES_URL', 'CPLES_XML', 'CPLE_AppDefined'
5 , 'CPLE_AssertionFailed', 'CPLE_FileIO', 'CPLE_IllegalArg', 'CPLE_NoWriteAccess'
6 , 'CPLE_None', 'CPLE_NotSupported', 'CPLE_OpenFailed', 'CPLE_OutOfMemory', 'CPLE
7 _UserInterrupt', 'CXT_Attribute', 'CXT_Comment', 'CXT_Element', 'CXT_Literal', '
8 CXT_Text', 'DCAP_CREATE', 'DCAP_CREATECOPY', 'DMD_CREATIONDATATYPES', 'DMD_CREAT
9 IONOPTIONLIST', 'DMD_EXTENSION', 'DMD_HELPTOPIC', 'DMD_LONGNAME', 'DMD_MIMETYPE'
10 , 'GA_ReadOnly', 'GA_Update', 'GCI_AlphaBand', 'GCI_BlackBand', 'GCI_BlueBand',
11 'GCI_CyanBand', 'GCI_GrayIndex', 'GCI_GreenBand', 'GCI_HueBand', 'GCI_LightnessB
12 and', 'GCI_MagentaBand', 'GCI_PaletteIndex', 'GCI_RedBand', 'GCI_SaturationBand'
13 , 'GCI_Undefined', 'GCI_YellowBand', 'GDT_Byte', 'GDT_CFloat32', 'GDT_CFloat64',
14 'GDT_CInt16', 'GDT_CInt32', 'GDT_Float32', 'GDT_Float64', 'GDT_Int16', 'GDT_Int
15 32', 'GDT_TypeCount', 'GDT_UInt16', 'GDT_UInt32', 'GDT_Unknown', 'GF_Read', 'GF_
16 Write', 'GPI_CMYK', 'GPI_Gray', 'GPI_HLS', 'GPI_RGB', 'GRA_Bilinear', 'GRA_Cubic
17 ', 'GRA_CubicSpline', 'GRA_NearestNeighbour', '__builtins__', '__doc__', '__file
18 __', '__name__']
19 >>>
那些GDT开头的就是数值数据类型。
1 >>> band.GetNoDataValue()
2 65535.0
3 >>> band.GetMaximum()
4 >>> band.GetMinimum()
5 >>> band.ComputeRasterMinMax()
6 (1.0, 255.0)
7 >>>
