使用VRT裁剪栅格影像

代码同样来自《Python地理数据处理》，这里做一个学习记录。书籍原外文名“Geoprocessing with Python”，在曼宁出版社可以搜索到，可以线上免费预览全文，下载源码和配套数据。附一个链接：Geoprocessing with Python

VRT（GDAL虚拟栅格格式），允许使用XML文件来定义数据集，本身不存储数据，只描述了如何从其他文件中获得数据。我自己的理解觉得有点像数据库里视图与数据表中视图的概念。

1.实现代码和输出结果

# -*- coding: utf-8 -*-
"""
Created on Fri Mar 25 13:27:15 2022

@author: Administrator

使用VRT裁剪影像

#列表或元组前面加星号作用是将列表解开成两个独立的参数，传入函数
"""

import os
from osgeo import gdal

os.chdir(r'E:\pycode\DATA\osgeopy-data\Landsat\Washington')
tmp_ds = gdal.Open('nat_color.tif')
tmp_gt = tmp_ds.GetGeoTransform()
inv_gt = gdal.InvGeoTransform(tmp_gt)
if gdal.VersionInfo()[0] == '1':
    if inv_gt[0] ==1:
        inv_gt = inv_gt[1]
    else:
        raise RuntimeError('Inverse geotransform failed')
elif inv_gt is None:
    raise RuntimeError('Inverse geotransform failed')
    
vashon_ul = (532000, 5262600)
vashon_lr = (548500, 5241500)
ulx, uly = map(int, gdal.ApplyGeoTransform(inv_gt, *vashon_ul))
lrx, lry = map(int, gdal.ApplyGeoTransform(inv_gt, *vashon_lr))
rows = lry - uly
columns = lrx - ulx
gt = list(tmp_gt)
gt[0] += gt[1] * ulx#新的原点x,y值
gt[3] += gt[5] * uly

ds = gdal.GetDriverByName('vrt').Create('vashon.vrt', columns, rows, 3)
ds.SetProjection(tmp_ds.GetProjection())
ds.SetGeoTransform(gt)

xml = '''
<SimpleSource>
    <SourceFilename relativeToVRT="1"> {fn} </SourceFilename>
    <SourceBand> {band} </SourceBand>
    <SrcRect xOff="{xoff}" yOff="{yoff}" xSize="{cols}" ySize="{rows}" />
    <DstRect xOff="0" yOff="0" xSize="{cols}" ySize="{rows}" />
</SimpleSource>
'''

data = {'fn': 'nat_color.tif', 
        'band': 1,
        'xoff': ulx,
        'yoff': uly,
        'cols': columns,
        'rows': rows
        }

meta = {'source_0': xml.format(**data)}#对字典解包，分别传入参数
ds.GetRasterBand(1).SetMetadata(meta, 'vrt_sources')

data['band'] = 2
meta = {'source_0': xml.format(**data)}
ds.GetRasterBand(2).SetMetadata(meta, 'vrt_sources')

data['band'] = 3
meta = {'source_0': xml.format(**data)}
ds.GetRasterBand(3).SetMetadata(meta, 'vrt_sources')

del ds, tmp_ds

生成的vashon.vrt文件，是用XML语法编写的，直接在浏览器中打开可以查看细节。

 QGIS软件支持打开VRT，上述文件打开效果：

 2.实验过程中犯了一个小的拼写错误，记录一下。

之前尝试在QGIS中打开生成的washon.vrt，一直报错，后来发现是自己敲代码时relativeToVRT敲错了字母。。。导致生成文件这个属性一直等于0，即relativeToVRT="0"了，找不到引用的源数据，.vrt就啥也不是。

此外，需注意，在能正常打开的情况，源数据文件nat_color.tif与虚拟栅格格式文件vashon.vrt必须是在同一目录下的。