矢栅处理(3)：Python读写GeoTIFF文件

发表于 分类于 阅读次数： Waline： 本文字数： 3.6k 阅读时长 ≈ 3 分钟

先看一下tif文件的数值类型和压缩方式。地理栅格数据通常会以GeoTIFF(.tif; .tiff)、HDF(.h5; .hdf5)、NETCDF(.nc)等多种文件格式(后两者气象数据居多，长时间序列数据存储)，但实际中空间处理、统计及可视化最常用的还是geotiff数据，因此先小结并对比一下：(1) gdal模块和rasterio两个模块读写tif的对比；(2)geotif文件的压缩格式

  • 关于h5文件和nc文件的读写后续总结

GeoTIFF文件

文件简介

GEO+TIF，如果 TIFF 文件本身没有地理配准信息，GDAL 将依次检查并使用扩展名为.aux.xml文件 . MapInfo .tab 文件以及tfw、.tifw/.tiffw 或 .wld 的 ESRI世界文件。

tiff和png等本身即是光栅文件，geotiff文件是在tiff文件中嵌入地理参考信息/tags得到的，包含两部分

  • 像元值，纯数值矩阵numpy ndarray
    • 矩阵的维度——行列数+波段数
    • 矩阵的数值及numpy数据类型dtype, 表示为，np.uint8; np.bool_; np.complex64
      • int8, int16, int32, int64; uint16, uint16, uint32, uint64;(无符号即大于等于0)
      • float16, float32, float64
  • 地理参考/空间位置信息
    • 空间范围Extent，矩形框四个角的范围
    • 坐标系：地理坐标系及投影坐标系，如果深究坐标参考系的信息较为丰富
    • 分辨率，或行列数、空值 Nodata value、图层个数等
  • tiff等光栅文件均可视为点文件，以该点为中心的格网空间范围内的值皆=该点的值
  • 数据读取和处理等中间过程可通过numpy数组完成，涉及空间运算时把位置信息加入生成进行处理即可

文件读写

简要版，更多的参数设置查看文档：rasterio; gdal1, gdal2. 仅以单波段图像为例进行读写分析，多波段可以自行拓展

读写后的tif信息往往不完全一致，可以通过读取两幅影像比较数值矩阵是否相等,下述三种方法写出文件的结果数值一致

1
(data1 == data2).all() #True则表示相同

rasterio(推荐)

1
2
3
4
5
6
7
8
9
10
11
12
13
import rasterio
# 读
with rasterio.open(file) as src_dataset:
    profiles = src_dataset.profile
    band1 = src_dataset.read(1)
print(profiles)
print(band1) #ndarray
# 在此更新字典信息
new_type = rasterio.int32
profiles.update(dtype=new_type,compress='DEFLATE')
# 写: band1为一个numpy array
with rasterio.open(out, 'w', **profiles) as dst:
        dst.write(band1.astype(new_type), 1)

gdal

更为底层，前者实际上是基于gdal完成的读写

gdal写出的文件会比rasterio更大，暂时还不知道原因

读，获取各项信息，可以通过dir(ds)查看其全部属性方法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
ds = gdal.Open(tiffile)
band = ds.GetRasterBand(1)
im_data = band.ReadAsArray()
print(im_data.dtype)
# 各项信息
print("文件的源信息",ds.GetMetadata())
im_width = ds.RasterXSize #栅格矩阵的列数
im_height = ds.RasterYSize #栅格矩阵的行数
im_bands = ds.RasterCount #波段数

im_geotrans = ds.GetGeoTransform() #仿射矩阵，左上角像素的大地坐标和像素分辨率
im_proj = ds.GetProjection() #地图投影信息，字符串表示
proj = osr.SpatialReference(wkt=im_proj)
nodata = band.GetNoDataValue()

写GDALDataset *Create(const char *pszName,nXSize,nYSize,nBands,eType,...)

GDALDataType, 如gdal.GDT_Byte

enum GDALDataType 对应numpy类型 备注
GDT_Byte np.uint8
GDT_UInt16, GDT_Int16 np.uint16, np.int16
GDT_UInt32, GDT_Int32 np.uint32, np.int32
GDT_UInt64, GDT_Int64 np.uint64, np.int64 GDAL >= 3.5
GDT_Float32， GDT_Float64 np.float32, np.float64
GDT_CInt16，GDT_CInt32 Complex Int16
GDT_CFloat32，GDT_CFloat64 Complex Float32

写出方法1: Create

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def write_geotif(outpath, im_data, im_geotrans,im_proj, nodata=-1):
    # 仅适用于单个波段
    # 数据类型
    GDT_Dict = {
        'uint8':gdal.GDT_Byte,
        'int8':gdal.GDT_Byte,
        'uint16':gdal.GDT_UInt16,
        'int16':gdal.GDT_Int16,
        'uint32':gdal.GDT_UInt32,
        'int32':gdal.GDT_Int32,
        'float32':gdal.GDT_Float32}
    gdaltype = GDT_Dict[im_data.dtype.name]
    im_bands, (im_height, im_width) = 1, im_data.shape

    # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(outpath, im_width, im_height, im_bands, gdaltype, options=["COMPRESS=LZW"])
    dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数
    dataset.SetProjection(im_proj) #写入投影
    # 写入数据
    dataset.GetRasterBand(1).WriteArray(im_data)
    dataset.GetRasterBand(1).SetNoDataValue(nodata)
    # 写入内存
    del dataset
outpath = r"....tif"
write_geotif(outpath, im_data, im_geotrans,im_proj, nodata=nodata)

写出方法2: CreateCopy

1
2
3
4
5
6
7
8
9
10
11

def write_geotif2(outpath, im_data, src_ds):
    # 创建文件
    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.CreateCopy(outpath, src_ds, strict=0, options=["COMPRESS=LZW"])
    # 写入数据
    dst_ds.GetRasterBand(1).WriteArray(im_data)
    # 写入内存
    del dst_ds
outpath = r"...tif"
write_geotif2(outpath, im_data,ds)

matlab

由于文件本质上是纯数值文件，Matlab处理tiff文件也十分便捷api，但鉴于目前不太熟悉，后续填充

参考

  1. https://rasterio.readthedocs.io/en/latest/intro.html
  2. https://gdal.org/tutorials/index.html
  3. https://gdal.org/index.html
  4. https://gdal.org/drivers/raster/gtiff.html
  5. 其他参考