『GDAL』读写栅格影像数据

2020-09-10

GDAL

字数统计: 1k | 阅读时长≈ 4 分钟

使用 GDAL读写tiff遥感影像数据

常用方法

打开已有的栅格数据集

fn = 'aster.img'
ds = gdal.Open(fn, GA_ReadOnly)
if ds is None:
   print 'Could not open ' + fn
    sys.exit(1)

读取栅格数据集的x方向像素数，y方向像素数，和波段数
```
cols = ds.RasterXSize
rows = ds.RasterYSize
bands = ds.RasterCount
```
注意后面没有括号，因为他们是属性(properties)不是方法(methods) 读取地理坐标参考信息(georeference info)
GeoTransform是一个list，存储着栅格数据集的地理坐标信息
```
adfGeoTransform[0] /* top left x 左上角x坐标*/
adfGeoTransform[1] /* w--e pixel resolution 东西方向上的像素分辨率*/
adfGeoTransform[2] /* rotation, 0 if image is "north up" 如果北边朝上，地图的旋转角度*/
adfGeoTransform[3] /* top left y 左上角y坐标*/
adfGeoTransform[4] /* rotation, 0 if image is "north up" 如果北边朝上，地图的旋转角度*/
adfGeoTransform[5] /* n-s pixel resolution 南北方向上的像素分辨率*/
```
注意栅格数据集的坐标一般都是以左上角为基准的。

下面的例子是从一个栅格数据集中取出Geotransform作为一个list，然后读取其中的数据
```
geotransform = ds.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
```
计算某一坐标对应像素的相对位置(pixel offset)，也就是该坐标与左上角的像素的相对位置，按像素数计算，计算公式如下：
```
xOffset = int((x – originX) / pixelWidth)
yOffset = int((y – originY) / pixelHeight)
```
读取某一像素点的值，需要分两步首先读取一个波段(band)：GetRasterBand()，其参数为波段的索引号然后用ReadAsArray(, , , )，读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。如果将矩阵大小设为1X1，就是读取一个像素了。但是这一方法只能将读出的数据放到矩阵中，就算只读取一个像素也是一样。例如：
```
band = ds.GetRasterBand(1)
data = band.ReadAsArray(xOffset, yOffset, 1, 1)
```
如果想一次读取一整张图，那么将offset都设定为0，size则设定为整个图幅的size，例如：
```
data = band.ReadAsArray(0, 0, cols, rows)
```
但是要注意，从data中读取某一像素的值，必须要用data[yoff, xoff]。注意不要搞反了。数学中的矩阵是[row,col]，而这里恰恰相反！这里面row对应y轴，col对应x轴。

读图像文件

from osgeo import gdal
# 读图像文件
def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
    im_proj = dataset.GetProjection()  # 地图投影信息
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height).astype(np.float)  # 将数据写成数组，对应栅格矩阵
    del dataset  # 关闭对象，文件dataset
    return im_proj, im_geotrans, im_data, im_height, im_width

完整程序

from osgeo import gdal
# 读图像文件
def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
    im_proj = dataset.GetProjection()  # 地图投影信息
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height).astype(np.float)  # 将数据写成数组，对应栅格矩阵
    del dataset  # 关闭对象，文件dataset
    return im_proj, im_geotrans, im_data, im_height, im_width

def write_img(filename, im_proj, im_geotrans, im_data):
    # gdal数据类型包括
    # gdal.GDT_Byte,
    # gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
    # gdal.GDT_Float32, gdal.GDT_Float64

    # 判断栅格数据的数据类型
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    # 判读数组维数
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape

    # 创建文件
    driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影

    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset
if __name__=='__main__':
    proj, geotrans, values, row, column = read_img(输入数据)  # 读数据
     write_img(r'输出地址', proj, geotrans, 输出影像名称)#写数据

参考

打赏

版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！