通过上一节介绍的方法，可以访问遥感影像的描述性信息，可以概括地知道影像的获取时间、处理时间、空间分辨率、影像大小等一些信息。 但是为了对遥感影像进行处理，需要进一步访问遥感影像中的数据，即影像中像元的灰度值。
GDAL提供了下面两个函数来访问影像的数值。
ReadRaster()读取图像数据(以二进制的形式)ReadAsArray()读取图像数据(以数组的形式)
In [9]:
from osgeo import gdal
dataset = gdal.Open("/data/gdata/lu75c.tif")
help(dataset.ReadRaster)
Help on method ReadRaster in module osgeo.gdal: ReadRaster(xoff=0, yoff=0, xsize=None, ysize=None, buf_xsize=None, buf_ysize=None, buf_type=None, band_list=None, buf_pixel_space=None, buf_line_space=None, buf_band_space=None, resample_alg=0, callback=None, callback_data=None, buf_obj=None) method of osgeo.gdal.Dataset instance
上面的 help() 函数查看gdal方法。
继续查看 ReadAsArray() 方法。
In [10]:
help(dataset.ReadAsArray)
Help on method ReadAsArray in module osgeo.gdal:
ReadAsArray(xoff=0, yoff=0, xsize=None, ysize=None, buf_obj=None, buf_xsize=None, buf_ysize=None, buf_type=None, resample_alg=0, callback=None, callback_data=None, interleave='band', band_list=None) method of osgeo.gdal.Dataset instance
Reading a chunk of a GDAL band into a numpy array. The optional (buf_xsize,buf_ysize,buf_type)
parameters should generally not be specified if buf_obj is specified. The array is returned
这是两个非常重要的函数，它们可以直接读取图像的数据，从而对栅格数据进行分析。
可以看到两个函数有许多的参数。解释一下：
xoff,yoff：指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位）。xsize,ysize： 指定要读取部分图像的矩形的长和宽（以像元为单位）。buf_xsize，buf_ysize：可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高， GDAL 将帮你缩放到这个大小。buf_type：可以对读出的数据的类型进行转换（比如原图数据类型是short，你要把它们缩小成byte）。band_list：适应多波段的情况。可以指定要读取的波段。
这里简单看一下如何获取GeoTIFF文件中的数据。
In [11]:
import array
from numpy import *
dataset.ReadAsArray(2500,2500,3,3)
array([[12, 12, 12],
[12, 12, 12],
[12, 12, 12]],dtype=int16)
dataset.ReadRaster(2500,2500,3,3)
Out[11]:
bytearray(b'\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00\x0c\x00')
我们就把图像左上角位于(3340, 3780)，宽高都为10个像元的数据读取出来了。
这两个函数返回的结果不一样，其中ReadAsArray()读出的是numpy的数组，类型为int16；
而ReadData()读出的是二进制型，具体的解释见后面。
ReadAsArray()返回的结果。关于类型的更多解释，见下一节的表[tab:gdalconst]。
In [12]:
from osgeo import gdal
# from osgeo.gdalconst import *
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
band = dataset.GetRasterBand(1)
band.ReadAsArray(100,100,5,5,10,10)
Out[12]:
array([[236, 236, 237, 237, 237, 237, 237, 237, 227, 227],
[236, 236, 237, 237, 237, 237, 237, 237, 227, 227],
[235, 235, 232, 232, 233, 233, 234, 234, 225, 225],
[235, 235, 232, 232, 233, 233, 234, 234, 225, 225],
[242, 242, 235, 235, 232, 232, 233, 233, 224, 224],
[242, 242, 235, 235, 232, 232, 233, 233, 224, 224],
[254, 254, 244, 244, 238, 238, 237, 237, 229, 229],
[254, 254, 244, 244, 238, 238, 237, 237, 229, 229],
[246, 246, 246, 246, 248, 248, 250, 250, 235, 235],
[246, 246, 246, 246, 248, 248, 250, 250, 235, 235]], dtype=uint8)
ReadAsArray的原型
help(band.ReadAsArray)分析help(band.ReadAsArray)函数的几个参数的意义。 前两个100是取值窗口的左上角在实际数据中所处象元的(x, y)位置。 后两个是取值窗口覆盖的区域大小， 再后面是取值窗口取出数组进行缩放后数组的大小。 这里需要注意的是这里的buffer大小是根据参数自动分配的，可以不指定， 如果不指定，则和第3，4个参数一致。经过5，6两个参数的设置，可以进行缩放。 假如取值窗口大小是20 × 20，取出后数组就可以人工设置大小。 让它成为10 × 10的数组，或者是40 × 40的数组。 如果设置成20 × 40的数组则取出的数组对于真实图像来说有了变形。
band.ReadAsArray(100,100,10,10)通过几个例子可以看到，读取的4个size参数的作用就是把硬盘上指定区域 （前四个参数定义）的数据按比例缩放到用户指定区域（后两个定义）内， 必要的时候进行缩放。 这里需要注意的是缩的情况（缩的时候是取几个周围点的平均值） 如果是调色板数据就可能引发问题。
栅格数据范围的处理
现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办？
In [13]:
band.XSize
Out[13]:
1500
In [14]:
band.YSize
Out[14]:
900
In [15]:
band.ReadAsArray(95,100,5,5)
Out[15]:
array([[246, 242, 237, 233, 234],
[247, 242, 237, 234, 234],
[244, 244, 240, 239, 242],
[241, 249, 245, 248, 254],
[244, 251, 252, 251, 247]], dtype=uint8)
ERROR 5: Access window out of range in RasterIO(). Requested
(95,100) of size 5x5 on raster of 100x100.可以看到，出错了。 获取数据的时候不能越界， 调用的时候要判断越界没。 还好用Python的numpy模块可以很方便地扩展矩阵。
In [16]:
from osgeo import gdal
import time
dataset = gdal.Open("/data/gdata/lu75c.tif")
band = dataset.GetRasterBand(1)
width, height = dataset.RasterXSize, dataset.RasterYSize
bw, bh = 128, 128
bxsize = width/bw
bysize = height/bh
band.ReadAsArray(0,0,width,height)
start = time.time()
band.ReadAsArray(0,0,width,height)
print (time.time()-start)
start2 = time.time()
for i in range(int(bysize)):
for j in range(int(bxsize)):
band.ReadAsArray(bw*j,bh*i,bw,bh)
print (time.time()-start2)
0.025071382522583008 0.09862303733825684
最后一段的循环的两个for位置调换一下。
运行一下得到结果：
运行python ch03_test_read_tif_time.py，得到以下结果：
0.03625917434692383
0.06595683097839355
0.06917715072631836上面第9行，目的是为了先读取一遍数据，不然，第一次的结果会略大一些，可能会超过第二次，会让人以为把图像分块读取比不分块读取的效率要高。
至于出现这种情况的原因，在于影像的存储是有不同方式的。
In [17]:
from osgeo import gdal
import numpy
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
band2 = dataset.GetRasterBand(2)
band3 = dataset.GetRasterBand(3)
cols = 100
rows = 100
data2 = band2.ReadAsArray(0, 0, cols,rows).astype(numpy.float16)
data3 = band3.ReadAsArray(0, 0, cols,rows).astype(numpy.float16)
mask = numpy.greater(data3 + data2, 0)
ndvi = numpy.choose(mask, (-99, (data3 - data2) / (data3 + data2)))
/tmp/ipykernel_733/2271967309.py:11: RuntimeWarning: invalid value encountered in divide ndvi = numpy.choose(mask, (-99, (data3 - data2) / (data3 + data2)))
In [18]:
ndvi
Out[18]:
array([[0.2766, 0.318 , 0.4285, ..., 0.3057, 0.275 , 0.6 ],
[0.1951, 0.2208, 0.3135, ..., 0.2683, 0.3538, 1. ],
[0.2142, 0.2195, 0.2896, ..., 0.4182, 0.758 , 1. ],
...,
[0.6313, 0.6113, 0.6553, ..., 0.3333, 0.359 , 0.36 ],
[0.4119, 0.4119, 0.422 , ..., 0.3784, 0.3784, 0.3418],
[0.3684, 0.322 , 0.309 , ..., 0.4084, 0.4 , 0.3684]],
dtype=float16)
In [19]:
from osgeo import gdalconst
dir(gdalconst)
Out[19]:
['CE_Debug', 'CE_Failure', 'CE_Fatal', 'CE_None', 'CE_Warning', 'CPLES_BackslashQuotable', 'CPLES_CSV', 'CPLES_SQL', 'CPLES_SQLI', 'CPLES_URL', 'CPLES_XML', 'CPLES_XML_BUT_QUOTES', 'CPLE_AWSAccessDenied', 'CPLE_AWSBucketNotFound', 'CPLE_AWSInvalidCredentials', 'CPLE_AWSObjectNotFound', 'CPLE_AWSSignatureDoesNotMatch', 'CPLE_AppDefined', 'CPLE_AssertionFailed', 'CPLE_FileIO', 'CPLE_HttpResponse', 'CPLE_IllegalArg', 'CPLE_NoWriteAccess', 'CPLE_None', 'CPLE_NotSupported', 'CPLE_ObjectNull', 'CPLE_OpenFailed', 'CPLE_OutOfMemory', 'CPLE_UserInterrupt', 'CXT_Attribute', 'CXT_Comment', 'CXT_Element', 'CXT_Literal', 'CXT_Text', 'DCAP_COORDINATE_EPOCH', 'DCAP_CREATE', 'DCAP_CREATECOPY', 'DCAP_CREATECOPY_MULTIDIMENSIONAL', 'DCAP_CREATE_FIELD', 'DCAP_CREATE_LAYER', 'DCAP_CREATE_MULTIDIMENSIONAL', 'DCAP_CURVE_GEOMETRIES', 'DCAP_DEFAULT_FIELDS', 'DCAP_DELETE_FIELD', 'DCAP_DELETE_LAYER', 'DCAP_FEATURE_STYLES', 'DCAP_FEATURE_STYLES_READ', 'DCAP_FEATURE_STYLES_WRITE', 'DCAP_FIELD_DOMAINS', 'DCAP_FLUSHCACHE_CONSISTENT_STATE', 'DCAP_GNM', 'DCAP_MEASURED_GEOMETRIES', 'DCAP_MULTIDIM_RASTER', 'DCAP_MULTIPLE_VECTOR_LAYERS', 'DCAP_NONSPATIAL', 'DCAP_NOTNULL_FIELDS', 'DCAP_NOTNULL_GEOMFIELDS', 'DCAP_OPEN', 'DCAP_RASTER', 'DCAP_RELATIONSHIPS', 'DCAP_RENAME_LAYERS', 'DCAP_REORDER_FIELDS', 'DCAP_SUBCREATECOPY', 'DCAP_UNIQUE_FIELDS', 'DCAP_VECTOR', 'DCAP_VIRTUALIO', 'DCAP_Z_GEOMETRIES', 'DIM_TYPE_HORIZONTAL_X', 'DIM_TYPE_HORIZONTAL_Y', 'DIM_TYPE_PARAMETRIC', 'DIM_TYPE_TEMPORAL', 'DIM_TYPE_VERTICAL', 'DMD_ALTER_FIELD_DEFN_FLAGS', 'DMD_ALTER_GEOM_FIELD_DEFN_FLAGS', 'DMD_CONNECTION_PREFIX', 'DMD_CREATIONDATATYPES', 'DMD_CREATIONFIELDDATASUBTYPES', 'DMD_CREATIONFIELDDATATYPES', 'DMD_CREATIONOPTIONLIST', 'DMD_CREATION_FIELD_DEFN_FLAGS', 'DMD_CREATION_FIELD_DOMAIN_TYPES', 'DMD_EXTENSION', 'DMD_EXTENSIONS', 'DMD_GEOMETRY_FLAGS', 'DMD_HELPTOPIC', 'DMD_ILLEGAL_FIELD_NAMES', 'DMD_LONGNAME', 'DMD_MIMETYPE', 'DMD_MULTIDIM_ARRAY_CREATIONOPTIONLIST', 'DMD_MULTIDIM_ARRAY_OPENOPTIONLIST', 'DMD_MULTIDIM_ATTRIBUTE_CREATIONOPTIONLIST', 'DMD_MULTIDIM_DATASET_CREATIONOPTIONLIST', 'DMD_MULTIDIM_DIMENSION_CREATIONOPTIONLIST', 'DMD_MULTIDIM_GROUP_CREATIONOPTIONLIST', 'DMD_NUMERIC_FIELD_WIDTH_INCLUDES_DECIMAL_SEPARATOR', 'DMD_NUMERIC_FIELD_WIDTH_INCLUDES_SIGN', 'DMD_OPENOPTIONLIST', 'DMD_SUBDATASETS', 'DMD_SUPPORTED_SQL_DIALECTS', 'GARIO_COMPLETE', 'GARIO_ERROR', 'GARIO_PENDING', 'GARIO_UPDATE', 'GA_ReadOnly', 'GA_Update', 'GCI_AlphaBand', 'GCI_BlackBand', 'GCI_BlueBand', 'GCI_CyanBand', 'GCI_GrayIndex', 'GCI_GreenBand', 'GCI_HueBand', 'GCI_LightnessBand', 'GCI_MagentaBand', 'GCI_PaletteIndex', 'GCI_RedBand', 'GCI_SaturationBand', 'GCI_Undefined', 'GCI_YCbCr_CbBand', 'GCI_YCbCr_CrBand', 'GCI_YCbCr_YBand', 'GCI_YellowBand', 'GDAL_DATA_COVERAGE_STATUS_DATA', 'GDAL_DATA_COVERAGE_STATUS_EMPTY', 'GDAL_DATA_COVERAGE_STATUS_UNIMPLEMENTED', 'GDAL_DCAP_CREATE_RELATIONSHIP', 'GDAL_DCAP_DELETE_RELATIONSHIP', 'GDAL_DCAP_UPDATE_RELATIONSHIP', 'GDAL_DMD_RELATIONSHIP_FLAGS', 'GDAL_DMD_RELATIONSHIP_RELATED_TABLE_TYPES', 'GDT_Byte', 'GDT_CFloat32', 'GDT_CFloat64', 'GDT_CInt16', 'GDT_CInt32', 'GDT_Float32', 'GDT_Float64', 'GDT_Int16', 'GDT_Int32', 'GDT_Int64', 'GDT_Int8', 'GDT_TypeCount', 'GDT_UInt16', 'GDT_UInt32', 'GDT_UInt64', 'GDT_Unknown', 'GDsCAddRelationship', 'GDsCDeleteRelationship', 'GDsCUpdateRelationship', 'GFT_Integer', 'GFT_Real', 'GFT_String', 'GFU_Alpha', 'GFU_AlphaMax', 'GFU_AlphaMin', 'GFU_Blue', 'GFU_BlueMax', 'GFU_BlueMin', 'GFU_Generic', 'GFU_Green', 'GFU_GreenMax', 'GFU_GreenMin', 'GFU_Max', 'GFU_MaxCount', 'GFU_Min', 'GFU_MinMax', 'GFU_Name', 'GFU_PixelCount', 'GFU_Red', 'GFU_RedMax', 'GFU_RedMin', 'GF_Read', 'GF_Write', 'GMF_ALL_VALID', 'GMF_ALPHA', 'GMF_NODATA', 'GMF_PER_DATASET', 'GPI_CMYK', 'GPI_Gray', 'GPI_HLS', 'GPI_RGB', 'GRA_Average', 'GRA_Bilinear', 'GRA_Cubic', 'GRA_CubicSpline', 'GRA_Lanczos', 'GRA_Max', 'GRA_Med', 'GRA_Min', 'GRA_Mode', 'GRA_NearestNeighbour', 'GRA_Q1', 'GRA_Q3', 'GRA_RMS', 'GRA_Sum', 'GRC_MANY_TO_MANY', 'GRC_MANY_TO_ONE', 'GRC_ONE_TO_MANY', 'GRC_ONE_TO_ONE', 'GRIORA_Average', 'GRIORA_Bilinear', 'GRIORA_Cubic', 'GRIORA_CubicSpline', 'GRIORA_Gauss', 'GRIORA_Lanczos', 'GRIORA_Mode', 'GRIORA_NearestNeighbour', 'GRIORA_RMS', 'GRTT_ATHEMATIC', 'GRTT_THEMATIC', 'GRT_AGGREGATION', 'GRT_ASSOCIATION', 'GRT_COMPOSITE', 'GTO_BIT', 'GTO_BSQ', 'GTO_TIP', 'OF_ALL', 'OF_GNM', 'OF_MULTIDIM_RASTER', 'OF_RASTER', 'OF_READONLY', 'OF_SHARED', 'OF_UPDATE', 'OF_VECTOR', 'OF_VERBOSE_ERROR', '_SwigNonDynamicMeta', '__builtin__', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_gdalconst', '_swig_add_metaclass', '_swig_python_version_info', '_swig_repr', '_swig_setattr_nondynamic_class_variable', '_swig_setattr_nondynamic_instance_variable']
那些GDT开头的就是数值数据类型。
想要查看图像中某一波段的数据类型，只需要打印这一波段的DataType属性即可。
In [20]:
from osgeo import gdal
dataset = gdal.Open("/data/gdata/geotiff_file.tif")
band = dataset.GetRasterBand(1)
band.DataType
Out[20]:
1
返回结果为整型。原来1表示的是 gdalconst.GDT_Byte 。注意这里的类型是与numpy中的类型对应的。
下面我们来看一个gdalconst与整型的对应值。
- 未知或未指定类型
gdalconst.GDT_Unknown0 - 8位无符整型
gdalconst.GDT_Byte1 - 16位无符整型
gdalconst.GDT_UInt162 - 16位整型
gdalconst.GDT_Int163 - 32位无符整型
gdalconst.GDT_UInt324 - 32位整型值
gdalconst.GDT_Int325 - 32位浮点型
gdalconst.GDT_Float326 - 64位浮点型
gdalconst.GDT_Float647 - 16位复数整型
gdalconst.GDT_CInt168 - 32位复数整型
gdalconst.GDT_CInt329 - 32位复数浮点型
gdalconst.GDT_CFloat3210 - 64位复数浮点型
gdalconst.GDT_CFloat6411