from osgeo import gdal
dataset = gdal.Open("/gdata/lu75c.tif")
help(dataset.ReadRaster)
上面的 help() 函数查看gdal方法:退出查看页面在终端输入“q”。继续查看 ReadAsArray() 方法。
help(dataset.ReadAsArray)
同样需要键入“q”退出。
这是两个非常重要的函数,它们可以直接读取图像的数据, 从而对栅格数据进行分析。可以看到两个函数的帮助中有许多的参数。 解释一下:
xoff,yoff
:指定想要读取的部分原点位置在整张图像中距离全图原点的位置(以像元为单位)。xsize,ysize
: 指定要读取部分图像的矩形的长和宽(以像元为单位)。buf_xsize,buf_ysize
:可以在读取出一部分图像后进行缩放。那么就用这两个参数来定义缩放后图像最终的宽和高,gdal将帮你缩放到这个大小。buf_type
:可以对读出的数据的类型进行转换(比如原图数据类型是short,你要把它们缩小成byte)。band_list
:适应多波段的情况。可以指定要读取的波段。
这里简单看一下如何获取GeoTIFF文件中的数据。
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)
我们就把图像左上角位于(3340, 3780),宽高都为10个像元的数据读取出来了。
这两个函数返回的结果不一样,其中ReadAsArray()
读出的是numpy的数组,类型为int16;
而ReadData()
读出的是二进制型,具体的解释见后面。
ReadAsArray()
返回的结果。关于类型的更多解释,见下一节的表[tab:gdalconst]。
from osgeo import gdal
from gdalconst import *
dataset = gdal.Open("/gdata/geotiff_file.tif")
band = dataset.GetRasterBand(1)
band.ReadAsArray(100,100,5,5,10,10)
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参数的作用就是把硬盘上指定区域 (前四个参数定义)的数据按比例缩放到用户指定区域(后两个定义)内, 必要的时候进行缩放。 这里需要注意的是缩的情况(缩的时候是取几个周围点的平均值) 如果是调色板数据就可能引发问题。
栅格数据范围的处理
现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办?
band.XSize
band.YSize
band.ReadAsArray(95,100,5,5)
ERROR 5: Access window out of range in RasterIO(). Requested
(95,100) of size 5x5 on raster of 100x100.
可以看到,出错了。 获取数据的时候不能越界, 调用的时候要判断越界没。 还好用Python的numpy模块可以很方便地扩展矩阵。
from osgeo import gdal
import time
dataset = gdal.Open("/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)
最后一段的循环的两个for位置调换一下。
运行一下得到结果:
运行python ch03_test_read_tif_time.py,得到以下结果:
0.03625917434692383
0.06595683097839355
0.06917715072631836
上面第9行,目的是为了先读取一遍数据,不然,第一次的结果会略大一些,可能会超过第二次,会让人以为把图像分块读取比不分块读取的效率要高。
至于出现这种情况的原因,在于影像的存储是有不同方式的。
from osgeo import gdal
import numpy
dataset = gdal.Open("/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)))
ndvi
from osgeo import gdalconst
dir(gdalconst)[:6] + ['... ...'] + dir(gdalconst)[-3:]
那些GDT开头的就是数值数据类型。
想要查看图像中某一波段的数据类型,只需要打印这一波段的DataType属性即可。
from osgeo import gdal
dataset = gdal.Open("/gdata/geotiff_file.tif")
band = dataset.GetRasterBand(1)
band.DataType
返回结果为整型。原来1表示的是 gdalconst.GDT_Byte
。注意这里的类型是与numpy中的类型对应的。
下面我们来看一个gdalconst与整型的对应值。
- 未知或未指定类型
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