7. GDAL python教程(6)——更多栅格数据处理函数

栅格数据的投影

首先要知道输入投影和输出投影的WKT(Well Known Text),可以通过GetProjection()读到,可以用SpatialReference对象创建。

用下面的语句新建栅格数据集并重新投影,投影结果输出到新数据集之中:

gdal.CreateAndReprojectImage( <source_dataset>, <output_filename>, src_wkt=<source_wkt>,       dst_wkt=<output_wkt>, dst_driver=<Driver>, eResampleAlg = <GDALResampleAlg>)

下面关注一下性能问题吧

对比两种方法的执行效率:逐个像素处理,使用内建函数处理

内建的函数效率最高,比逐个像素循环要快很多很多。

7.1. 滑动窗口

7.1.1. 按坐标逐个循环

data = inBand.ReadAsArray(0, 0, cols, rows).astype(Numeric.Int)
outData = Numeric.zeros((rows, cols), Numeric.Float)
for i in range(1, rows-1): # skipping first & last
   for j in range(1, cols-1):
     outData[i,j] = (data[i-1,j-1] + data[i-1,j] + data[i-1,j+1] + data[i,j-1] + data[i,j] + data[i,j+1] + data[i+1,j-1] + data[i+1,j] + data[i+1,j+1]) / 9.0

7.1.2. 数组切片

data = inBand.ReadAsArray(0, 0, cols, rows).astype(Numeric.Int))
outData = Numeric.zeros((rows, cols), Numeric.Int)
outData[1:rows-1,1:cols-1] = (data[0:rows-2,0:cols-2] + data[0:rows-2,1:cols-1] + data[0:rows-2,2:cols] + data[1:rows-1, 0:cols-2] + data[1:rows-1,1:cols-1] + data[1:rows-1,2:cols] + data[2:rows,0:cols-2] + data[2:rows,1:cols-1] + data[2:rows,2:cols]) / 9

用数组切片的方法相对而言效率较高。

7.2. 过滤器filter

或者叫滤波器吧。低通滤波Low-pass filter,用于平滑数据。高通滤波High-pass filter,用于边缘增强。上面滑动窗口的两段程序就是低通滤波。