1. 目的
(1) 掌握如何利用GDAL读取和写入栅格数据
(2) 学习了解GDAL的常用函数
2. 版本
2.1 2019年11月06日 Version 1
3. 任务简述
本示例中运用的是Landsat数据,任务目的是读取3景单一波段的Landsat影像,然后合并成一个堆栈影像。
4. 代码示例
# 1. 基本模块导入
import os
from osgeo import gdal
# 2. 路径处理和基本变量定义
os.chdir(r'F:\01_Data')
band1_fn = 'p047r027_7t20000730_z10_nn10.tif'
band2_fn = 'p047r027_7t20000730_z10_nn20.tif'
band3_fn = 'p047r027_7t20000730_z10_nn30.tif'
# 3. 基于band1的属性创建三波段GeoTIFF
in_ds = gdal.Open(band1_fn)
print(in_ds)
in_band = in_ds.GetRasterBand(1) #打开波段1
print('*******************************')
print(in_band)
gtiff_driver = gdal.GetDriverByName('Gtiff') # Gtiff是待读取的数据的格式
# (1) 创建三波段Geotiff
out_ds = gtiff_driver.Create('Nat_Color2.tif',in_band.XSize,in_band.YSize,3,in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection()) #投影
out_ds.SetGeoTransform(in_ds.GetGeoTransform()) #地理坐标变换
# 4.从输入波段复制像素到输出波段3
in_data = in_band.ReadAsArray()
out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)
# 5.从数据集而不是波段复制像素
in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())
out_ds.GetRasterBand(1).WriteArray(gdal.Open(band3_fn).ReadAsArray())
# 6. 为每个波段计算统计值
out_ds.FlushCache()
for i in range(1,4):
out_ds.GetRasterBand(i).ComputeStatistics(False)
out_ds.BuildOverviews('average',[2,4,8,16,32]) #创建视图或金字塔图层
del out_ds
print('Finish')
5. 参考资料
[1] Python地理数据处理
[2] Python的gdal.GetDriverByName实例