Python地理数据处理 十三：栅格数据处理（一）

1.数据重采样

  重采样是指根据一类象元的信息内插出另一类象元信息的过程。在遥感中，重采样是从高分辨率遥感影像中提取出低分辨率影像的过程。常用的重采样方法有最邻近内插法、双线性内插法和三次卷积法内插。
  ReadAsArray函数可以重采样读取的数据，并且指定输出缓冲区大小或传递一个已有的缓冲区数组。
函数格式：

band.ReadAsArray([xoff], [yoff], [win_xsize], [win_ysize], [buf_xsize], [buf_ysize], [buf_obj])

1. xoff是开始阅读的专栏，默认值为0。
2. yoff是开始阅读的行，默认值为0。
3. win_xsize是要读取的列数，默认为全部读取。
4. win_ysize是要读取的行数，默认为全部读取。
5. buf_xsize是输出数组中的列数，默认值为使用win_xsize的值。
6. buf_ysize是输出数组中的行数，默认值为使用win_ysize的值。
7. buf_obj是一个NumPy数组，用于将数据放入其中，而不是创建一个新数组。如果需要，数据将被重新采样以适合此数组，对应的值也将转换为该数组的数据类型。

  具有较大尺寸的数组将重采样为较小的像素，更小的尺寸的数组，将采用最近邻插值法重采样为更大的像素

重采样为更小的像素：

  需要提供一个比原有数据维数更大的数组，使得重复使用像素值来填充目标数组：

import os
from osgeo import gdal

os.chdir(r'D:\geodata\Landsat\Washington')

in_ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
in_band = in_ds.GetRasterBand(1)

# 计算输出行列数
# 输入数翻倍，因为我将像素大小减半
out_rows = in_band.YSize * 2
out_columns = in_band.XSize * 2

# 创建输出数据集
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('band1_resampled.tif',
    out_columns, out_rows)

# 编辑地理变换
# 像素变为原来的 1/4
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] /= 2
geotransform[5] /= 2
out_ds.SetGeoTransform(geotransform)

# 读取数据时，指定一个较大的缓冲
data = in_band.ReadAsArray(
    buf_xsize=out_columns, buf_ysize=out_rows)

# 将数据写入输出光栅
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)

# 为较大的图像构建合适的概视图
out_band.FlushCache()
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32, 64])

del out_ds

  图像展示：

  属性信息：

2. 字节序列

1.ReadRaster([xoff], [yoff], [xsize], [ysize], [buf_xsize], [buf_ysize], [buf_type], [band_list]，[buf_pixel_space], [buf_line_space], [buf_band_space])

1. xoff是列读取起点，默认值为0。
2. yoff是行读取起点，默认值为0。
3. xsize是读取的列数，默认为全部读取。
4. ysize是读取的行数，默认为全部读取。
5. buf_xsize是输出数组里的列数，默认值为使用 win_xsize 值。如果此值不同于 win_xsize，则数据将重采样。
6. buf_ysize是输出数组里的行数,默认值为使用 win_ysize 值。如果此值不同于 win_ysize，则数据将重采样。
7. buf_type 是返回序列的GDAL目标数据类型，默认值与源数据相同。
8. band_list是要读取的波段列表，默认读取所有波段。
9. buf_pixel_space是序列中像素之间的字节偏移，默认值为buf_type的大小。
10. buf_line_space是序列中行间的字节偏移，默认值是 buf_type 乘以 xsize 。
11. buf_band_space是序列中列间的字节偏移，默认值是 buf_line_space 乘以 ysize 。

import os
import numpy as np
from osgeo import gdal

data_dir = r'D:\geodata'

os.chdir(os.path.join(data_dir, 'Landsat', 'Washington'))
ds = gdal.Open('Landsat_color.tif')
data = ds.ReadRaster(1400, 6000, 2, 2, band_list=[1])
print(data)

# 取出第一个值，将从字节转换为数字
print(data[0])

# 尝试更改第一个像素的值。
# 输出失败，因为不能更改字节字符串
data[0] = 50

# 将字节字符串转换为字节数组
# 更改第一个值，输出成功
bytearray_data = bytearray(data)  # bytearray是字节数组
bytearray_data[0] = 50
print(bytearray_data[0])

b'\x1c\x1d\x1c\x1e'
28
50

  2. 将字节字符串转换为像素值的元组：

import struct
tuple_data = struct.unpack('B' * 4, data) # 指定4个字节
print(tuple_data)

(28, 29, 28, 30)

  3. 将元组转换为numpy数组：

numpy_data1 = np.array(tuple_data)
print(numpy_data1)

[28 29 28 30]

  4. 将字节字符串转换为numpy数组：

# 将字节字符串转换为numpy数组
# 重构一个numpy数组，使其具有2行2列，就像读入的原始数据一样
numpy_data2 = np.fromstring(data, np.int8)
reshaped_data = np.reshape(numpy_data2, (2,2))
print(reshaped_data)

[[28 29]
 [28 30]]

  5. 利用字节序列将图像重采样为更大的像素尺寸：

import os
import numpy as np
from osgeo import gdal

os.chdir(r'D:\geodata\Landsat\Washington')

# 打开输入栅格
in_ds = gdal.Open('Landsat_color.tif')

# 计算输出的行和列的数量
# 输入数字的一半，因为要使像素的两倍大
out_rows = int(in_ds.RasterYSize / 2)
out_columns = int(in_ds.RasterXSize / 2)
num_bands = in_ds.RasterCount

# 创建输出数据集
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('Landsat_color_resampled.tif',
    out_columns, out_rows, num_bands)

# 编辑地理变换，让像素尺寸变大
# 像素宽度、像素高度加倍
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] *= 2
geotransform[5] *= 2
out_ds.SetGeoTransform(geotransform)

# 利用较小的缓冲来读写数据
data = in_ds.ReadRaster(
    buf_xsize=out_columns, buf_ysize=out_rows)

# 将数据写入栅格中
out_ds.WriteRaster(0, 0, out_columns, out_rows, data)

# 为较小的图像创建合适数量的概视图
out_ds.FlushCache()
for i in range(num_bands):
    out_ds.GetRasterBand(i + 1).ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16])

del out_ds

结果显示：

3. 子数据集

  MODIS图像（中分辨率成像光谱仪）是一个分层数据格式（HDF）文件。如果数据包含子数据集，可用 GetSubDatasets() 函数得到他们的子数据集列表，再利用这个列表信息打开所需要的子数据集。

  数据集结构：

import os
import numpy as np
from osgeo import gdal

data_dir = r'D:\geodata'

# 从MODIS文件中获取子数据集的名称和描述
# 并打印，输出NDVI（归一化植被指数）
os.chdir(os.path.join(data_dir, 'Modis'))
ds = gdal.Open('MYD13Q1.A2014313.h20v11.005.2014330092746.hdf')
subdatasets = ds.GetSubDatasets() # GetSubDatasets函数返回一个元组列表
print('Number of subdatasets: {}'.format(len(subdatasets)))
for sd in subdatasets:
    print('Name: {0}\nDescription:{1}\n'.format(*sd))


# 打开Modis文件中的第一个子数据集：[0][0]
# 第五个数据集为: [4][0]
ndvi_ds = gdal.Open(subdatasets[0][0])

# 通过打印尺寸来确保正常工作
print('Dataset dimensions: {} {}'.format(ndvi_ds.RasterXSize, ndvi_ds.RasterYSize))

# 读取数据之前先获取波段
ndvi_band = ndvi_ds.GetRasterBand(1)
print('Band dimensions: {} {}'.format(ndvi_band.XSize, ndvi_band.YSize))

结果显示：