VRT—GDAL虚拟格式

司机简称

VRT

Driver built-in by default

This driver is built-in by default

介绍

VRT驱动程序是GDAL的一个格式驱动程序,它允许虚拟GDAL数据集由其他GDAL数据集组成,并具有重新定位、可能应用的算法以及更改或添加的各种元数据。数据集的VRT描述可以保存为XML格式,通常给出扩展名.VRT。

VRT格式还可以描述 扭曲VRT平底锅

An example of a simple .vrt file referring to a 512x512 dataset with one band loaded from utm.tif might look like this:

<VRTDataset rasterXSize="512" rasterYSize="512">
    <GeoTransform>440720.0, 60.0, 0.0, 3751320.0, 0.0, -60.0</GeoTransform>
    <VRTRasterBand dataType="Byte" band="1">
        <ColorInterp>Gray</ColorInterp>
        <SimpleSource>
        <SourceFilename relativeToVRT="1">utm.tif</SourceFilename>
        <SourceBand>1</SourceBand>
        <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
        <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
        </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

VRT文件的许多方面都是 栅格数据模型 为了理解不同元素的语义,应该对其进行审查。

VRT文件可以通过转换为VRT格式生成。然后可以编辑结果文件以修改映射、添加元数据或其他目的。VRT文件也可以通过各种方式以编程方式生成。

本教程将介绍.vrt文件格式(适合用户编辑.vrt文件),以及如何为开发人员以编程方式创建和操作.vrt文件。

.vrt格式

A XML schema of the GDAL VRT format is available.

存储在磁盘上的虚拟文件以XML格式保存,包含以下元素。

VRTDataset :这是整个GDAL数据集的根元素。它必须具有rasterXSize和rasterYSize属性,用像素描述数据集的宽度和高度。它可能具有值为VRTWarpedDataset的子类属性 (扭曲VRT )或VRTPansharpenedDataset (平底锅 ). 它可能有SRS、GeoTransform、GCPList、Metadata、MaskBand和VRTRasterBand子元素。

<VRTDataset rasterXSize="512" rasterYSize="512">

VRTDataset

VRTDataset允许的子元素有:

  • SRS: This element contains the spatial reference system (coordinate system) in OGC WKT format. Note that this must be appropriately escaped for XML, so items like quotes will have the ampersand escape sequences substituted. As well as WKT, valid input to the OGRSpatialReference::SetFromUserInput() method (such as well known GEOGCS names, and PROJ.4 format) is also allowed in the SRS element.

<SRS dataAxisToSRSAxisMapping="1,2">PROJCS[&quot;NAD27 / UTM zone 11N&quot;,GEOGCS[&quot;NAD27&quot;,DATUM[&quot;North_American_Datum_1927&quot;,SPHEROID[&quot;Clarke 1866&quot;,6378206.4,294.9786982139006,AUTHORITY[&quot;EPSG&quot;,&quot;7008&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6267&quot;]],PRIMEM[&quot;Greenwich&quot;,0],UNIT[&quot;degree&quot;,0.0174532925199433],AUTHORITY[&quot;EPSG&quot;,&quot;4267&quot;]],PROJECTION[&quot;Transverse_Mercator&quot;],PARAMETER[&quot;latitude_of_origin&quot;,0],PARAMETER[&quot;central_meridian&quot;,-117],PARAMETER[&quot;scale_factor&quot;,0.9996],PARAMETER[&quot;false_easting&quot;,500000],PARAMETER[&quot;false_northing&quot;,0],UNIT[&quot;metre&quot;,1,AUTHORITY[&quot;EPSG&quot;,&quot;9001&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;26711&quot;]]</SRS>

这个 数据轴到rsaxisMapping 自GDAL 3.0以来,允许属性描述CRS定义中指示的轴与GeoTransform或GCP元数据的轴之间的关系。属性的值是以逗号分隔的整数列表。此列表的元素数必须是CRS的轴数。值从1开始。如果m表示此属性的数组值,则m [0] 是CRS第一个轴的数据轴号。如果属性丢失,则隐含OAMS_TRADITIONAL_GIS_ORDER data axis to CRS axis mapping策略。

  • GeoTransform :此元素包含数据集的六值仿射地理变换,即像素/线坐标与地理参考坐标之间的映射。

<GeoTransform>440720.0,  60,  0.0,  3751320.0,  0.0, -60.0</GeoTransform>
  • GCPList :此元素包含数据集的地面控制点列表、像素/线坐标和地理参考坐标之间的映射。“投影”属性应包含地理参考坐标的SRS,格式与SRS元素相同。dataAxisToSRSAxisMapping属性与SRS元素中的属性相同。

<GCPList Projection="EPSG:4326">
    <GCP Id="1" Info="a" Pixel="0.5" Line="0.5" X="0.0" Y="0.0" Z="0.0" />
    <GCP Id="2" Info="b" Pixel="13.5" Line="23.5" X="1.0" Y="2.0" Z="0.0" />
</GCPList>
  • 元数据 :此元素包含与整个VRTDataset或VRTRasterBand关联的元数据名称/值对的列表。它有<MDI>(元数据项)子元素,这些子元素具有“key”属性和作为元素数据的值。元数据元素可以重复多次,在这种情况下,它必须附带一个“domain”属性以指示元数据域的名称。

<Metadata>
  <MDI key="md_key">Metadata value</MDI>
</Metadata>
  • MaskBand :此元素表示数据集上所有带之间共享的掩码带(请参阅RFC 15中的GMFU PERU数据集)。它必须包含单个VRTRasterBand子元素,即掩码带本身的描述。

<MaskBand>
  <VRTRasterBand dataType="Byte">
    <SimpleSource>
      <SourceFilename relativeToVRT="1">utm.tif</SourceFilename>
      <SourceBand>mask,1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
      <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
    </SimpleSource>
  </VRTRasterBand>
</MaskBand>
  • OverviewList: (GDAL >= 3.2.0, not valid for VRTPansharpenedDataset) This elements contains a list of overview factors, separated by space, to create "virtual overviews". For example 2 4. It can be used so that bands of the VRT datasets declare overviews. This only makes sense to use if the sources added in those bands have themselves overviews compatible with the declared factor. It is generally not necessary to use this mechanism, since downsampling pixel requests on a VRT dataset/band are able to use overviews of the sources, even when the VRT bands do not declare them. One situation where explicit overviews are needed at the VRT level is the warping of a VRT to a lower resolution. This element can also be used with an existing VRT dataset by running GDALDataset::BuildOverviews() or gdaladdo with the VRT_VIRTUAL_OVERVIEWS configuration option set to YES. Virtual overviews have the least priority compared to the Overview element at the VRTRasterBand level, or to materialized .vrt.ovr files.

  • VRTRasterBand :这表示数据集的一个频带。

VRTRasterBand

VRTRasterBand的属性包括:

  • 数据库类型 (可选):与此波段关联的像素数据的类型(使用名称Byte、UInt16、Int16、UInt32、Int32、Float32、Float64、CInt16、CInt32、CFloat32或CFloat64)。如果未指定,则默认为1

  • band (可选):此元素表示的带区编号(基于1)。

  • 方块大小 (可选,GDAL>=3.3):块宽度。如果未指定,则默认为栅格宽度和128之间的最小值。

  • 块状化 (可选,GDAL>=3.3):块高度。如果未指定,则默认为栅格高度的最小值。

此元素可能有Metadata、ColorInterp、NoDataValue、HideNoDataValue、ColorTable、gdalMasterAttributeTable、Description和MaskBand子元素以及各种源元素,如SimpleSource、ComplexSource等。栅格带可能有许多“源”,指示实际栅格数据应从何处获取,以及如何将其映射到栅格带像素空间。

VRTRasterBand允许的子元素有:

  • ColorInterp :此元素的数据应该是颜色解释类型的名称。灰色、调色板、红色、绿色、蓝色、Alpha、色调、饱和度、亮度、青色、洋红、黄色、黑色或未知色之一。

<ColorInterp>Gray</ColorInterp>:
  • NoDataValue: If the input datasets to be composed have a nodata value for this raster band, set this element's value to that nodata value for it to be reflected in the VRT. This must not be confused with the NODATA element of a VRTComplexSource element.

<NoDataValue>-100.0</NoDataValue>
  • HideNoDataValue :如果此值为1,则不会报告nodata值。本质上,调用方在读取nodata像素时不会意识到它。从中复制/转换的任何数据集都没有nodata值。当您要为数据集指定固定的背景值时,这非常有用。背景将是NoDataValue元素指定的值。当此元素不存在时,默认值为0。

<HideNoDataValue>1</HideNoDataValue>
  • ColorTable :此元素是定义颜色表中的项的一组项元素的父元素。目前只支持RGBA颜色表,c1为红色,c2为绿色,c3为蓝色,c4为alpha。条目是有序的,并假定从颜色表条目0开始。

<ColorTable>
  <Entry c1="0" c2="0" c3="0" c4="255"/>
  <Entry c1="145" c2="78" c3="224" c4="255"/>
</ColorTable>
  • GDALRasterAttributeTable: (GDAL >=2.3) This element is parent to a set of FieldDefn elements defining the columns of a raster attribute table, followed by a set of Row elements defining the values of the columns of each row.

<GDALRasterAttributeTable>
  <FieldDefn index="0">
    <Name>Value</Name>
    <Type>0</Type>
    <Usage>0</Usage>
  </FieldDefn>
  <FieldDefn index="1">
    <Name>Red</Name>
    <Type>0</Type>
    <Usage>6</Usage>
  </FieldDefn>
  <FieldDefn index="2">
    <Name>Green</Name>
    <Type>0</Type>
    <Usage>7</Usage>
  </FieldDefn>
  <FieldDefn index="3">
    <Name>Blue</Name>
    <Type>0</Type>
    <Usage>8</Usage>
  </FieldDefn>
  <Row index="0">
    <F>-500</F>
    <F>127</F>
    <F>40</F>
    <F>65</F>
  </Row>
  <Row index="1">
    <F>-400</F>
    <F>154</F>
    <F>168</F>
    <F>118</F>
  </Row>
</GDALRasterAttributeTable>
  • 描述 :此元素包含栅格标注栏的可选说明作为其文本值。

<Description>Crop Classification Layer</Description>
  • UnitType :此可选元素包含高程标注栏数据的垂直单位。米是“m”,英尺是“ft”。默认假设为米。

<UnitType>ft</UnitType>
  • 抵消 :此可选元素包含从栅格带上缩放的像素值计算“实际”像素值时应应用的偏移量。默认值为0.0。

<Offset>0.0</Offset>
  • 规模 :此可选元素包含当从栅格带上的缩放像素值计算“真实”像素值时应应用的比例。默认值为1.0。

<Scale>0.0</Scale>
  • 概述 :此可选元素描述带区的一个概述级别。它应该有一个子SourceFilename和SourceBand元素。SourceFilename可能具有relativeToVRT布尔属性。可以使用多个元素来描述多个概述。

<Overview>
  <SourceFilename relativeToVRT="1">yellowstone_2.1.ntf.r2</SourceFilename>
  <SourceBand>1</SourceBand>
</Overview>
  • CategoryNames :此可选元素包含类别子元素列表,其中包含已分类栅格标注栏的类别名称。

<CategoryNames>
  <Category>Missing</Category>
  <Category>Non-Crop</Category>
  <Category>Wheat</Category>
  <Category>Corn</Category>
  <Category>Soybeans</Category>
</CategoryNames>
  • SimpleSource: The SimpleSource indicates that raster data should be read from a separate dataset, indicating the dataset, and band to be read from, and how the data should map into this band's raster space.

  • AveragedSource :AveragedSource从SimpleSource派生,并共享相同的属性,只是当目标矩形的大小与源矩形的大小不同时,它使用平均重采样而不是SimpleSource中的最近邻算法。注意:可以使用更通用的机制来指定重采样算法。请参阅上面有关“重采样”属性的段落。

  • ComplexSource: The ComplexSource is derived from the SimpleSource (so it shares the SourceFilename, SourceBand, SrcRect and DstRect elements), but it provides support to rescale and offset the range of the source values. Certain regions of the source can be masked by specifying the NODATA value, or starting with GDAL 3.3, with the <UseMaskBand>true</UseMaskBand> element.

  • KernelFilteredSource: The KernelFilteredSource is a pixel source derived from the Simple Source (so it shares the SourceFilename, SourceBand, SrcRect and DstRect elements, but it also passes the data through a simple filtering kernel specified with the Kernel element.

  • MaskBand :此元素表示特定于它所包含的VRTRasterBand的掩码带。它必须包含单个VRTRasterBand子元素,即掩码带本身的描述。

来源

SimpleSource

SimpleSource可以有SourceFilename、SourceBand、srcret和DstRect子元素。srrect元素将指示应读取指定源文件上的哪个矩形,DstRect元素指示应如何将源数据的矩形映射到VRTRasterBands空间。

SourceFilename上的relativeToVRT属性指示文件名是应解释为相对于.vrt文件(值为1)还是不应解释为相对于.vrt文件(值为0)。默认值为0。

Some characteristics of the source band can be specified in the optional SourceProperties element to enable the VRT driver to defer the opening of the source dataset until it really needs to read data from it. This is particularly useful when building VRTs with a big number of source datasets. The needed parameters are the raster dimensions, the size of the blocks and the data type. If the SourceProperties tag is not present, the source dataset will be opened at the same time as the VRT itself.

备注

Starting with GDAL 3.4, the SourceProperties element is no longer necessary for deferred opening of the source datasets.

SourceBand子元素的内容可以引用掩码带。例如,mask,1表示源的第一个频带的掩码频带。

<SimpleSource>
  <SourceFilename relativeToVRT="1">utm.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="512" RasterYSize="512" DataType="Byte" BlockXSize="128" BlockYSize="128"/>
  <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
  <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
</SimpleSource>

可以添加OpenOptions子元素来指定打开源数据集时要应用的打开选项。它有<OOI>(打开选项项)子元素,这些子元素有一个“key”属性和作为元素数据的值。

<SimpleSource>
  <SourceFilename relativeToVRT="1">utm.tif</SourceFilename>
  <OpenOptions>
      <OOI key="OVERVIEW_LEVEL">0</OOI>
  </OpenOptions>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="256" RasterYSize="256" DataType="Byte" BlockXSize="128" BlockYSize="128"/>
  <SrcRect xOff="0" yOff="0" xSize="256" ySize="256"/>
  <DstRect xOff="0" yOff="0" xSize="256" ySize="256"/>
</SimpleSource>

可以在SimpleSource或ComplexSource元素上指定重采样属性,以指定当目标矩形的大小与源矩形的大小不同时使用的重采样算法。该属性允许的值有:nearest、bilinear、cubic、cubicspline、lanczos、average和mode。

<SimpleSource resampling="cubic">
  <SourceFilename relativeToVRT="1">utm.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="256" RasterYSize="256" DataType="Byte" BlockXSize="128" BlockYSize="128"/>
  <SrcRect xOff="0" yOff="0" xSize="256" ySize="256"/>
  <DstRect xOff="0" yOff="0" xSize="128" ySize="128"/>
</SimpleSource>
ComplexSource

除了线性缩放之外,还可以通过指定Exponent、SrcMin、SrcMax、DstMin和DstMax元素来使用使用幂函数的非线性缩放。如果未指定SrcMin和SrcMax,则根据源最小值和最大值计算它们(这可能需要分析整个源数据集)。指数必须为正。(这5个值可以通过gdalu translate的-exponent和-scale选项进行设置。)

ComplexSource支持添加自定义查找表以将源值转换为目标值。可以使用以下形式指定LUT:

<LUT>[src value 1]:[dest value 1],[src value 2]:[dest value 2],...</LUT>

中间值使用对应范围的边界目标值之间的线性插值计算。

ComplexSource支持从具有颜色表的源栅格带中提取颜色组件。ColorTableComponent值是要提取的颜色组件的索引:1表示红色带,2表示绿色带,3表示蓝色带,或4表示alpha带。

转换源值时,按以下顺序执行操作:

  • 掩蔽,如果NODATA元素被设置,或者从gdal3.3开始,如果UseMaskBand被设置为true并且源频带具有掩蔽频带。请注意,这只是二进制掩蔽,因此如果掩蔽带实际上是具有非0或非255值的alpha带,则不会进行alpha混合。

  • 颜色表扩展

  • 对于线性缩放,应用缩放比率,然后缩放偏移

  • 对于非线性缩放,应用(DstMax DstMin)*pow((SrcValue SrcMin)/(SrcMax SrcMin),指数)+DstMin

  • 查表

<ComplexSource>
  <SourceFilename relativeToVRT="1">utm.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <ScaleOffset>0</ScaleOffset>
  <ScaleRatio>1</ScaleRatio>
  <ColorTableComponent>1</ColorTableComponent>
  <LUT>0:0,2345.12:64,56789.5:128,2364753.02:255</LUT>
  <NODATA>0</NODATA>  <!-- if the mask is a mask or alpha band, use <UseMaskBand>true</UseMaskBand> -->
  <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
  <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
</ComplexSource>

非线性缩放:

<ComplexSource>
  <SourceFilename relativeToVRT="1">16bit.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <Exponent>0.75</Exponent>
  <SrcMin>0</SrcMin>
  <SrcMax>65535</SrcMax>
  <DstMin>0</DstMin>
  <DstMax>255</DstMax>
  <SrcRect xOff="0" yOff="0" xSize="512" ySize="512"/>
  <DstRect xOff="0" yOff="0" xSize="512" ySize="512"/>
</ComplexSource>
KernelFilteredSource

内核元素应该有两个子元素,Size和Coefs以及可选的规范化布尔属性(默认值为false=0)。大小必须始终是奇数,并且Coefs必须具有由空格分隔的size*size条目。目前,核不适用于欠采样或过采样数据。

<KernelFilteredSource>
  <SourceFilename>/debian/home/warmerda/openev/utm.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <Kernel normalized="1">
    <Size>3</Size>
    <Coefs>0.11111111 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111 0.11111111</Coefs>
  </Kernel>
</KernelFilteredSource>

从GDAL 2.3开始,也可以使用可分离的内核。在这种情况下,Coefs条目的数量应该与大小相对应。Coefs指定了一个一维内核,该内核沿着每个轴依次应用,从而使执行速度更快。许多常用的图像处理滤波器是可分离的。例如,高斯模糊:

<KernelFilteredSource>
  <SourceFilename>/debian/home/warmerda/openev/utm.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <Kernel normalized="1">
    <Size>13</Size>
    <Coefs>0.01111 0.04394 0.13534 0.32465 0.60653 0.8825 1.0 0.8825 0.60653 0.32465 0.13534 0.04394 0.01111</Coefs>
  </Kernel>
</KernelFilteredSource>

概览

当处理涉及下采样的RasterIO()请求时,GDAL可以有效地利用组成波段的源中可用的概述。但在一般情况下,VRT乐队本身不会暴露概述。

除非(从最高优先级到较低优先级):

  • 这个 概述 元素存在于VRTRasterBand元素中。见上文。

  • 或外部。vrt.ovr概述已构建

  • (从gdal3.2开始)显式虚拟概述,如果 OverviewList 元素在VRTDataset元素中声明(见上文)。这些虚拟的概述将被外部隐藏。vrt.ovr公司以后可能构建的概述。

  • (从gdal2.1开始)隐式虚拟概述,如果VRTRasterBand由具有概述的单个SimpleSource或ComplexSource组成。这些虚拟的概述将被外部隐藏。vrt.ovr公司以后可能构建的概述。

原始文件的.vrt说明

到目前为止,我们已经描述了如何从GDAL支持的现有文件中派生新的虚拟数据集。但是,通常还需要使用原始二进制栅格文件,这些文件的数据常规布局已知,但不存在特定格式的驱动程序。这可以通过编写描述原始文件的.vrt文件来实现。

例如,下面的.vrt描述了一个原始栅格文件,其中包含名为<i>l2p3hhsso.img</i>的文件中的浮点复杂像素。图像数据从第一个字节开始(ImageOffset=0)。像素之间的字节偏移量为8(PixelOffset=8),即CFloat32的大小。从一行开始到下一行开始的字节偏移量为9376字节(line offset=9376),即宽度(1172)乘以像素(8)的大小。

<VRTDataset rasterXSize="1172" rasterYSize="1864">
    <VRTRasterBand dataType="CFloat32" band="1" subClass="VRTRawRasterBand">
        <SourceFilename relativetoVRT="1">l2p3hhsso.img</SourceFilename>
        <ImageOffset>0</ImageOffset>
        <PixelOffset>8</PixelOffset>
        <LineOffset>9376</LineOffset>
        <ByteOrder>MSB</ByteOrder>
    </VRTRasterBand>
</VRTDataset>

需要注意的是,VRTRawRasterBand有一个子类说明符“VRTRawRasterBand”。此外,VRTRawRasterBand包含许多以前未看到的元素,但没有“源”信息。vrtrawrasterband可能永远没有源(即SimpleSource),但除了前面描述的所有仍然受支持的常规“元数据”元素外,还应包含以下元素。

  • SourceFilename :包含此标注栏数据的原始文件的名称。relativeToVRT属性可用于指示源文件名是否与.vrt文件(1)相关(0)。

  • ImageOffset :此图像波段数据的第一个像素开始的偏移量(字节)。默认为零。

  • PixelOffset :从一个像素开始到同一行上下一个像素的偏移量(字节)。在压缩的单波段数据中,这将是 数据库类型 以字节为单位。

  • LineOffset :从一个数据扫描行开始到下一个数据扫描行的偏移量(字节)。在压缩的单波段数据中,这将是PixelOffset*rasterXSize。

  • ByteOrder :定义磁盘上数据的字节顺序。LSB(最低有效字节优先),如英特尔x86系统上的自然字节顺序;MSB(最高有效字节优先),如摩托罗拉或Sparc系统上的自然字节顺序。默认为本地机器订单。

其他一些注意事项:

  • 假设磁盘上的图像数据与频带的数据类型相同 数据库类型 在维特拉瓦斯特乐队。

  • 支持VRTRasterBand的所有非源属性,包括颜色表、元数据、nodata值和颜色解释。

  • VRTRawRasterBand支持栅格的就地更新,而基于源的VRTRasterBand始终是只读的。

  • OpenEV工具包含一个文件菜单选项,用于在GUI中输入描述原始栅格文件的参数并创建相应的.vrt文件。

  • 一个.vrt文件中的多个波段可以来自同一个原始文件。只需确保每个波段的ImageOffset、PixelOffset和LineOffset定义适合该特定波段的像素。

另一个例子,在本例中是400x300 RGB像素的交错图像。

<VRTDataset rasterXSize="400" rasterYSize="300">
<VRTRasterBand dataType="Byte" band="1" subClass="VRTRawRasterBand">
    <ColorInterp>Red</ColorInterp>
    <SourceFilename relativetoVRT="1">rgb.raw</SourceFilename>
    <ImageOffset>0</ImageOffset>
    <PixelOffset>3</PixelOffset>
    <LineOffset>1200</LineOffset>
</VRTRasterBand>
<VRTRasterBand dataType="Byte" band="2" subClass="VRTRawRasterBand">
    <ColorInterp>Green</ColorInterp>
    <SourceFilename relativetoVRT="1">rgb.raw</SourceFilename>
    <ImageOffset>1</ImageOffset>
    <PixelOffset>3</PixelOffset>
    <LineOffset>1200</LineOffset>
</VRTRasterBand>
<VRTRasterBand dataType="Byte" band="3" subClass="VRTRawRasterBand">
    <ColorInterp>Blue</ColorInterp>
    <SourceFilename relativetoVRT="1">rgb.raw</SourceFilename>
    <ImageOffset>2</ImageOffset>
    <PixelOffset>3</PixelOffset>
    <LineOffset>1200</LineOffset>
</VRTRasterBand>
</VRTDataset>

VRT数据集的创建

The VRT driver supports several methods of creating VRT datasets. The vrtdataset.h include file should be installed with the core GDAL include files, allowing direct access to the VRT classes. However, even without that most capabilities remain available through standard GDAL interfaces.

To create a VRT dataset that is a clone of an existing dataset use the GDALDriver::CreateCopy() method. For example to clone utm.tif into a wrk.vrt file in C++ the following could be used:

GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
GDALDataset *poSrcDS, *poVRTDS;

poSrcDS = (GDALDataset *) GDALOpenShared( "utm.tif", GA_ReadOnly );

poVRTDS = poDriver->CreateCopy( "wrk.vrt", poSrcDS, FALSE, NULL, NULL, NULL );

GDALClose((GDALDatasetH) poVRTDS);
GDALClose((GDALDatasetH) poSrcDS);

Note the use of GDALOpenShared() when opening the source dataset. It is advised to use GDALOpenShared() in this situation so that you are able to release the explicit reference to it before closing the VRT dataset itself. In other words, in the previous example, you could also invert the 2 last lines, whereas if you open the source dataset with GDALOpen(), you'd need to close the VRT dataset before closing the source dataset.

To obtain the resulting VRT XML of wrk.vrt without having to read the text from an actual file, you can modify the above code to open the new dataset with an empty filename and use the "xml:VRT" metadata domain.

// no filename
poVRTDS = poDriver->CreateCopy( "", poSrcDS, FALSE, NULL, NULL, NULL );

// obtain the actual XML text that a VRT file would contain
const char *xmlvrt = poVRTDS->GetMetadata("xml:VRT")[0];

To create a virtual copy of a dataset with some attributes added or changed such as metadata or coordinate system that are often hard to change on other formats, you might do the following. In this case, the virtual dataset is created "in memory" only by virtual of creating it with an empty filename, and then used as a modified source to pass to a GDALDriver::CreateCopy() written out in TIFF format.

poVRTDS = poDriver->CreateCopy( "", poSrcDS, FALSE, NULL, NULL, NULL );

poVRTDS->SetMetadataItem( "SourceAgency", "United States Geological Survey");
poVRTDS->SetMetadataItem( "SourceDate", "July 21, 2003" );

poVRTDS->GetRasterBand( 1 )->SetNoDataValue( -999.0 );

GDALDriver *poTIFFDriver = (GDALDriver *) GDALGetDriverByName( "GTiff" );
GDALDataset *poTiffDS;

poTiffDS = poTIFFDriver->CreateCopy( "wrk.tif", poVRTDS, FALSE, NULL, NULL, NULL );

GDALClose((GDALDatasetH) poTiffDS);

In the above example the nodata value is set as -999. You can set the HideNoDataValue element in the VRT dataset's band using GDALRasterBand::SetMetadataItem() on that band.

poVRTDS->GetRasterBand( 1 )->SetMetadataItem( "HideNoDataValue" , "1" );

In this example a virtual dataset is created with the GDALDriver::Create() method, and adding bands and sources programmatically, but still via the "generic" API. A special attribute of VRT datasets is that sources can be added to the VRTRasterBand (but not to VRTRawRasterBand) by passing the XML describing the source into GDALRasterBand::SetMetadataItem() on the special domain target "new_vrt_sources". The domain target "vrt_sources" may also be used, in which case any existing sources will be discarded before adding the new ones. In this example we construct a simple averaging filter source instead of using the simple source.

// construct XML for simple 3x3 average filter kernel source.
const char *pszFilterSourceXML  =
"<KernelFilteredSource>"
"  <SourceFilename>utm.tif</SourceFilename><SourceBand>1</SourceBand>"
"  <Kernel>"
"    <Size>3</Size>"
"    <Coefs>0.111 0.111 0.111 0.111 0.111 0.111 0.111 0.111 0.111</Coefs>"
"  </Kernel>"
"</KernelFilteredSource>";

// Create the virtual dataset.
poVRTDS = poDriver->Create( "", 512, 512, 1, GDT_Byte, NULL );
poVRTDS->GetRasterBand(1)->SetMetadataItem("source_0",pszFilterSourceXML,
                                            "new_vrt_sources");

A more general form of this that will produce a 3x3 average filtered clone of any input datasource might look like the following. In this case we deliberately set the filtered datasource as in the "vrt_sources" domain to override the SimpleSource created by the cpp:func:`GDALDriver::CreateCopy method. The fact that we used cpp:func:GDALDriver::CreateCopy ensures that all the other metadata, georeferencing and so forth is preserved from the source dataset ... the only thing we are changing is the data source for each band.

int   nBand;
GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
GDALDataset *poSrcDS, *poVRTDS;

poSrcDS = (GDALDataset *) GDALOpenShared( pszSourceFilename, GA_ReadOnly );

poVRTDS = poDriver->CreateCopy( "", poSrcDS, FALSE, NULL, NULL, NULL );

for( nBand = 1; nBand <= poVRTDS->GetRasterCount(); nBand++ )
{
    char szFilterSourceXML[10000];

    GDALRasterBand *poBand = poVRTDS->GetRasterBand( nBand );

    sprintf( szFilterSourceXML,
        "<KernelFilteredSource>"
        "  <SourceFilename>%s</SourceFilename><SourceBand>%d</SourceBand>"
        "  <Kernel>"
        "    <Size>3</Size>"
        "    <Coefs>0.111 0.111 0.111 0.111 0.111 0.111 0.111 0.111 0.111</Coefs>"
        "  </Kernel>"
        "</KernelFilteredSource>",
        pszSourceFilename, nBand );

    poBand->SetMetadataItem( "source_0", szFilterSourceXML, "vrt_sources" );
}

The VRTDataset class is one of the few dataset implementations that supports the GDALDataset::AddBand() method. The options passed to the GDALDataset::AddBand() method can be used to control the type of the band created (VRTRasterBand, VRTRawRasterBand, VRTDerivedRasterBand), and in the case of the VRTRawRasterBand to set its various parameters. For standard VRTRasterBand, sources should be specified with the above GDALRasterBand::SetMetadataItem() examples.

GDALDriver *poDriver = (GDALDriver *) GDALGetDriverByName( "VRT" );
GDALDataset *poVRTDS;

poVRTDS = poDriver->Create( "out.vrt", 512, 512, 0, GDT_Byte, NULL );
char** papszOptions = NULL;
papszOptions = CSLAddNameValue(papszOptions, "subclass", "VRTRawRasterBand"); // if not specified, default to VRTRasterBand
papszOptions = CSLAddNameValue(papszOptions, "SourceFilename", "src.tif"); // mandatory
papszOptions = CSLAddNameValue(papszOptions, "ImageOffset", "156"); // optional. default = 0
papszOptions = CSLAddNameValue(papszOptions, "PixelOffset", "2"); // optional. default = size of band type
papszOptions = CSLAddNameValue(papszOptions, "LineOffset", "1024"); // optional. default = size of band type * width
papszOptions = CSLAddNameValue(papszOptions, "ByteOrder", "LSB"); // optional. default = machine order
papszOptions = CSLAddNameValue(papszOptions, "relativeToVRT", "true"); // optional. default = false
poVRTDS->AddBand(GDT_Byte, papszOptions);
CSLDestroy(papszOptions);

delete poVRTDS;

使用派生的带(C/C++中的像素函数)

一种特殊的波段类型是“派生”波段,它从源波段派生其像素信息。对于这种类型的标注栏,还必须指定一个像素函数,该函数负责生成输出栅格。像素函数由应用程序创建,然后使用唯一键向GDAL注册。

Using derived bands you can create VRT datasets that manipulate bands on the fly without having to create new band files on disk. For example, you might want to generate a band using four source bands from a nine band input dataset (x0, x3, x4, and x8) and some constant y:

band_value = sqrt((x3*x3+x4*x4)/(x0*x8)) + y;

您可以编写pixel函数来计算这个值,然后用名为“MyFirstFunction”的GDAL注册它。然后,可以使用以下VRT XML来显示此派生频带:

<VRTDataset rasterXSize="1000" rasterYSize="1000">
    <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
        <Description>Magnitude</Description>
        <PixelFunctionType>MyFirstFunction</PixelFunctionType>
        <PixelFunctionArguments y="4" />
        <SimpleSource>
            <SourceFilename relativeToVRT="1">nine_band.dat</SourceFilename>
            <SourceBand>1</SourceBand>
            <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
            <DstRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
        </SimpleSource>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">nine_band.dat</SourceFilename>
            <SourceBand>4</SourceBand>
            <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
            <DstRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
        </SimpleSource>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">nine_band.dat</SourceFilename>
            <SourceBand>5</SourceBand>
            <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
            <DstRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
        </SimpleSource>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">nine_band.dat</SourceFilename>
            <SourceBand>9</SourceBand>
            <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
            <DstRect xOff="0" yOff="0" xSize="1000" ySize="1000"/>
        </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

备注

PixelFunctionArguments can only be used with C++ pixel functions in GDAL versions 3.4 and greater.

除了子类规范(VRTDerivedRasterBand)和PixelFunctionType值之外,还有一个新参数可以派上用场:SourceTransferType。通常,源栅格是使用导出频带的数据类型获得的。但是,有时您希望pixel函数能够访问比生成的数据类型更高分辨率的源数据。例如,您可能有一个“Float”类型的派生带,它接受“cfloatt32”或“cfloatt64”类型的单个源,并返回虚部。为此,请将SourceTransferType设置为“CFloat64”。否则,在调用像素函数之前,源将被转换为“Float”,虚部将丢失。

<VRTDataset rasterXSize="1000" rasterYSize="1000">
    <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
        <Description>Magnitude</Description>
        <PixelFunctionType>MyFirstFunction</PixelFunctionType>
        <SourceTransferType>CFloat64</SourceTransferType>
        ...

默认像素函数

GDAL provides a set of default pixel functions that can be used without writing new code:

PixelFunctionType

Number of input sources

PixelFunctionArguments

描述

cmul

2

multiply the first band for the complex conjugate of the second

complex

2

make a complex band merging two bands used as real and imag values

conj

1

computes the complex conjugate of a single raster band (just a copy if the input is non-complex)

dB

1

fact (optional)

perform conversion to dB of the abs of a single raster band (real or complex): 20. * log10( abs( x ) ). The optional fact parameter can be set to 10 to get the alternative formula: 10. * log10( abs( x ) )

dB2amp

1

perform scale conversion from logarithmic to linear (amplitude) (i.e. 10 ^ ( x / 20 ) ) of a single raster band (real only). Deprecated in GDAL v3.5. Please use the exp pixel function with base = 10. and fact = 0.05 i.e. 1./20

dB2pow

1

perform scale conversion from logarithmic to linear (power) (i.e. 10 ^ ( x / 10 ) ) of a single raster band (real only). Deprecated in GDAL v3.5. Please use the exp pixel function with base = 10. and fact = 0.1 i.e. 1./10

diff

2

computes the difference between 2 raster bands (b1 - b2)

div

2

divide one rasted band by another (b1 / b2)

exp

1

base (optional), fact (optional)

computes the exponential of each element in the input band x (of real values): e ^ x. The function also accepts two optional parameters: base and fact that allow to compute the generalized formula: base ^ ( fact * x ). Note: this function is the recommended one to perform conversion form logarithmic scale (dB): `` 10. ^ (x / 20.)``, in this case base = 10. and fact = 0.05 i.e. 1. / 20

imag

1

extract imaginary part from a single raster band (0 for non-complex)

intensity

1

computes the intensity Re( x * conj(x) ) of a single raster band (real or complex)

interpolate_exp

>= 2

t0, dt, t

interpolate a value at time (or position) t given input sources beginning at position t0 with spacing dt using exponential interpolation

interpolate_linear

>= 2

t0, dt, t

interpolate a value at time (or position) t given input sources beginning at t0 with spacing dt using linear interpolation

inv

1

k (optional)

inverse (1./x). If the optional k parameter is set then the result is multiplied by k (k / x)

log10

1

compute the logarithm (base 10) of the abs of a single raster band (real or complex): log10( abs( x ) )

mod

1

extract module from a single raster band (real or complex)

mul

>= 2

k (optional)

multiply 2 or more raster bands. If the optional k parameter is provided then the result is multiplied by the scalar k.

phase

1

extract phase from a single raster band [-PI,PI] (0 or PI for non-complex)

polar

2

amplitude_type (optional)

make a complex band using input bands for amplitude and phase values b1 * exp( j * b2 ). The optional (string) parameter amplitude_type can be AMPLITUDE (default) INTENSITY or dB. Note: if amplitude_type is set to INTENSITY then negative values are clipped to zero.

pow

1

power

raise a single raster band to a constant power, specified with argument power (real only)

real

1

extract real part from a single raster band (just a copy if the input is non-complex)

sqrt

1

perform the square root of a single raster band (real only)

sum

>= 2

k (optional)

sum 2 or more raster bands. If the optional k parameter is provided then it is added to each element of the result

replace_nodata

= 1

to (optional)

convert incoming NoData values to a new value, IEEE 754 nan by default

scale

= 1

perform scaling according to the offset and scale values of the raster band

写入像素函数

To register this function with GDAL (prior to accessing any VRT datasets with derived bands that use this function), an application calls GDALAddDerivedBandPixelFuncWithArgs() with a key and a GDALDerivedPixelFuncWithArgs:

static const char pszMetadata[] =
"<PixelFunctionArgumentsList>"
"   <Argument name='y' description='y' type='double' mandatory='1' />"
"   <Argument type='builtin' value='offset' />"
"   <Argument type='builtin' value='scale' />"
"   <Argument type='builtin' value='NoData' />"
"   <Argument name='customConstant' type='constant' value='42'>"
"</PixelFunctionArgumentsList>";
GDALAddDerivedBandPixelFuncWithArgs("MyFirstFunction", TestFunction, pszMetadata);

A good time to do this is at the beginning of an application when the GDAL drivers are registered. pszMetadata is optional and can be nullptr. It can be used to declare the function signature to the user and to request additional parameters aside from the ones from the Dataset.

A GDALDerivedPixelFuncWithArgs is defined with a signature similar to GDALRasterBand::IRasterIO():

CPLErr TestFunction(void **papoSources, int nSources, void *pData, int nBufXSize, int nBufYSize, GDALDataType eSrcType, GDALDataType eBufType, int nPixelSpace, int nLineSpace, CSLConstList papszArgs)
参数
  • papoSources -- A pointer to packed rasters; one per source. The datatype of all will be the same, specified in the eSrcType parameter.

  • nSources -- The number of source rasters.

  • pData -- The buffer into which the data should be read, or from which it should be written. This buffer must contain at least nBufXSize * nBufYSize words of type eBufType. It is organized in left to right, top to bottom pixel order. Spacing is controlled by the nPixelSpace and nLineSpace parameters.

  • nBufXSize -- The width of the buffer image into which the desired region is to be read, or from which it is to be written.

  • nBufYSize -- The height of the buffer image into which the desired region is to be read, or from which it is to be written.

  • eSrcType -- The type of the pixel values in the papoSources raster array.

  • eBufType -- The type of the pixel values that the pixel function must generate in the pData data buffer.

  • nPixelSpace -- The byte offset from the start of one pixel value in pData to the start of the next pixel value within a scanline. If defaulted (0) the size of the datatype eBufType is used.

  • nLineSpace -- The byte offset from the start of one scanline in pData to the start of the next.

  • papszArgs -- An optional string list of named function arguments (e.g. y=4)

It is also possible to register a GDALDerivedPixelFunc (which omits the final CSLConstList argument) using GDALAddDerivedBandPixelFunc().

以下是像素功能的实现:

#include "gdal.h"

CPLErr TestFunction(void **papoSources, int nSources, void *pData,
                    int nXSize, int nYSize,
                    GDALDataType eSrcType, GDALDataType eBufType,
                    int nPixelSpace, int nLineSpace,
                    CSLConstList papszArgs)
{
    int ii, iLine, iCol;
    double pix_val;
    double x0, x3, x4, x8;

    // ---- Init ----
    if (nSources != 4) return CE_Failure;

    const char *pszY = CSLFetchNameValue(papszArgs, "y");
    if (pszY == nullptr) return CE_Failure;

    double NoData = NAN;
    const char *pszNoData = CSLFetchNameValue(papszArgs, "NoData");
    if (pszNoData != nullptr)
    {
        NoData = std::strtod(pszNoData, &end);
        if (end == pszNoData) return CE_Failure; // Could not parse
    }

    char *end = nullptr;
    double y = std::strtod(pszY, &end);
    if (end == pszY) return CE_Failure; // Could not parse

    // ---- Set pixels ----
    for( iLine = 0; iLine < nYSize; iLine++ )
    {
        for( iCol = 0; iCol < nXSize; iCol++ )
        {
            ii = iLine * nXSize + iCol;
            /* Source raster pixels may be obtained with SRCVAL macro */
            x0 = SRCVAL(papoSources[0], eSrcType, ii);
            x3 = SRCVAL(papoSources[1], eSrcType, ii);
            x4 = SRCVAL(papoSources[2], eSrcType, ii);
            x8 = SRCVAL(papoSources[3], eSrcType, ii);

            if (x0 == NoData || x3 == NoData || x4 == NoData || x8 == NoData)
                pix_val = NAN;
            else
                pix_val = sqrt((x3*x3+x4*x4)/(x0*x8)) + y;

            GDALCopyWords(&pix_val, GDT_Float64, 0,
                        ((GByte *)pData) + nLineSpace * iLine + iCol * nPixelSpace,
                        eBufType, nPixelSpace, 1);
        }
    }

    // ---- Return success ----
    return CE_None;
}

使用派生带(在Python中使用像素函数)

Starting with GDAL 2.2, in addition to pixel functions written in C/C++ as documented in the 使用派生的带(C/C++中的像素函数) section, it is possible to use pixel functions written in Python. Both CPython and NumPy are requirements at run-time.

VRTRasterBand的子元素(其子类规范必须设置为VRTDerivedRasterBand)是:

  • PixelFunctionType (必需):必须设置为函数名,该函数名将在PixelFunctionCode元素中定义为内联Python模块,或以“module_name.function_name”的形式表示外部Python模块中的函数

  • PixelFunctionLanguage (必需):必须设置为Python。

  • PixelFunctionCode (如果PixelFunctionType的形式为“function_name”,则为必需,否则忽略)。Python模块的内嵌代码,它必须至少有一个名为PixelFunctionType的函数。

  • BufferRadius (可选,默认为0):相对于要满足的原始RasterIO()请求,在传递给像素函数的输入和输出缓冲区的左侧、右侧、底部和顶部获取的额外像素量。请注意,将忽略此缓冲区中输出缓冲区的值。

Python pixel函数的签名必须具有以下参数:

  • in_ar :输入NumPy数组列表(每个源一个NumPy数组)

  • out_ar :要填充的输出NumPy数组。数组以正确的维度初始化,并使用VRTRasterBand.dataType。

  • xoff :带区访问区域左上角的像素偏移量。通常不需要,除非处理取决于栅格中的像素位置。

  • yoff 到标注栏访问区域左上角的线偏移。一般不需要。

  • 大小 :带区访问区域的宽度。可与外形一起使用 [1] 以确定请求的水平重采样率。

  • 伊西泽 :带区访问区域的高度。可与外形一起使用 [0] 以确定请求的垂直重采样率。

  • raster_xsize :栅格带的总计。一般不需要。

  • raster_ysize :栅格带的总计。一般不需要。

  • buf_radius :添加到inu ar/outu ar的左、右、上和下的缓冲区半径(以像素为单位)。这是可选的buffer radius元素的值,可以对其进行设置,以便将原始像素请求扩展给定数量的像素。

  • gt :地理变换。6个双值数组。

  • 关键字参数 :在PixelFunctionArguments中定义了用户参数的字典

提供 out_ar 必须就地修改数组。将忽略pixel函数当前返回的任何值。

备注

如果你想填补 out_ar from another array, use the out_ar[:] = ... 语法。

实例

将源文件的值乘以系数1.5的VRT

<VRTDataset rasterXSize="20" rasterYSize="20">
    <SRS>EPSG:26711</SRS>
    <GeoTransform>440720,60,0,3751320,0,-60</GeoTransform>
    <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
        <PixelFunctionType>multiply</PixelFunctionType>
        <PixelFunctionLanguage>Python</PixelFunctionLanguage>
        <PixelFunctionArguments factor="1.5"/>
        <PixelFunctionCode><![CDATA[
            import numpy as np
            def multiply(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                            raster_ysize, buf_radius, gt, **kwargs):
                factor = float(kwargs['factor'])
                out_ar[:] = np.round_(np.clip(in_ar[0] * factor,0,255))
            ]]>
        </PixelFunctionCode>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">byte.tif</SourceFilename>
        </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

添加2个(或更多)栅格的VRT

<VRTDataset rasterXSize="20" rasterYSize="20">
    <SRS>EPSG:26711</SRS>
    <GeoTransform>440720,60,0,3751320,0,-60</GeoTransform>
    <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
        <PixelFunctionType>add</PixelFunctionType>
        <PixelFunctionLanguage>Python</PixelFunctionLanguage>
        <PixelFunctionCode><![CDATA[
            import numpy as np
            def add(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                            raster_ysize, buf_radius, gt, **kwargs):
                np.round_(np.clip(np.sum(in_ar, axis = 0, dtype = 'uint16'),0,255),
                        out = out_ar)
            ]]>
        </PixelFunctionCode>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">byte.tif</SourceFilename>
        </SimpleSource>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">byte2.tif</SourceFilename>
        </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

使用外部库强制山体阴影的VRT

<VRTDataset rasterXSize="121" rasterYSize="121">
    <SRS>EPSG:4326</SRS>
    <GeoTransform>-80.004166666666663,0.008333333333333,0,
    44.004166666666663,0,-0.008333333333333</GeoTransform>
    <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
        <ColorInterp>Gray</ColorInterp>
        <SimpleSource>
            <SourceFilename relativeToVRT="1">n43.dt0</SourceFilename>
        </SimpleSource>
        <PixelFunctionLanguage>Python</PixelFunctionLanguage>
        <PixelFunctionType>hillshading.hillshade</PixelFunctionType>
        <PixelFunctionArguments scale="111120" z_factor="30" />
        <BufferRadius>1</BufferRadius>
        <SourceTransferType>Int16</SourceTransferType>
    </VRTRasterBand>
</VRTDataset>

使用hillshading.py:

# Licence: MIT
# Copyright 2016, Even Rouault
import math

def hillshade_int(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                        raster_ysize, radius, gt, z, scale):
    ovr_scale_x = float(out_ar.shape[1] - 2 * radius) / xsize
    ovr_scale_y = float(out_ar.shape[0] - 2 * radius) / ysize
    ewres = gt[1] / ovr_scale_x
    nsres = gt[5] / ovr_scale_y
    inv_nsres = 1.0 / nsres
    inv_ewres = 1.0 / ewres

    az = 315
    alt = 45
    degreesToRadians = math.pi / 180

    sin_alt = math.sin(alt * degreesToRadians)
    azRadians = az * degreesToRadians
    z_scale_factor = z / (8 * scale)
    cos_alt_mul_z_scale_factor = \
            math.cos(alt * degreesToRadians) * z_scale_factor
    cos_az_mul_cos_alt_mul_z_scale_factor_mul_254 = \
                254 * math.cos(azRadians) * cos_alt_mul_z_scale_factor
    sin_az_mul_cos_alt_mul_z_scale_factor_mul_254 = \
                254 * math.sin(azRadians) * cos_alt_mul_z_scale_factor
    square_z_scale_factor = z_scale_factor * z_scale_factor
    sin_alt_mul_254 = 254.0 * sin_alt

    for j in range(radius, out_ar.shape[0]-radius):
        win_line = in_ar[0][j-radius:j+radius+1,:]
        for i in range(radius, out_ar.shape[1]-radius):
            win = win_line[:,i-radius:i+radius+1].tolist()
            x = inv_ewres * ((win[0][0] + win[1][0] + win[1][0] + win[2][0])-\
                            (win[0][2] + win[1][2] + win[1][2] + win[2][2]))
            y = inv_nsres * ((win[2][0] + win[2][1] + win[2][1] + win[2][2])-\
                            (win[0][0] + win[0][1] + win[0][1] + win[0][2]))
            xx_plus_yy = x * x + y * y
            cang_mul_254 = (sin_alt_mul_254 - \
                (y * cos_az_mul_cos_alt_mul_z_scale_factor_mul_254 - \
                    x * sin_az_mul_cos_alt_mul_z_scale_factor_mul_254)) / \
                math.sqrt(1 + square_z_scale_factor * xx_plus_yy)
            if cang_mul_254 < 0:
                out_ar[j,i] = 1
            else:
                out_ar[j,i] = 1 + round(cang_mul_254)

def hillshade(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
            raster_ysize, radius, gt, **kwargs):
    z = float(kwargs['z_factor'])
    scale= float(kwargs['scale'])
    hillshade_int(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                raster_ysize, radius, gt, z, scale)

Python模块路径

When importing modules from inline Python code or when relying on out-of-line code (PixelFunctionType of the form "module_name.function_name"), you need to make sure the modules are accessible through the python path. Note that contrary to the Python interactive interpreter, the current path is not automatically added when used from GDAL. So you may need to define the PYTHONPATH environment variable if you get ModuleNotFoundError exceptions.

安全影响

The ability to run Python code potentially opens the door to many potential vulnerabilities if the user of GDAL may process untrusted datasets. To avoid such issues, by default, execution of Python pixel function will be disabled. The execution policy can be controlled with the GDAL_VRT_ENABLE_PYTHON configuration option, which can accept 3 values:

  • 是:所有VRT脚本都被认为是可信的,当涉及像素操作时,它们的Python像素函数将运行。

  • 否:所有VRT脚本都被认为是不可信的,并且不会运行任何PythonPixelFunction。

  • TRUSTED_MODULES (default setting): all VRT scripts with inline Python code in their PixelFunctionCode elements will be considered untrusted and will not be run. VRT scripts that use a PixelFunctionType of the form "module_name.function_name" will be considered as trusted, only if "module_name" is allowed in the GDAL_VRT_TRUSTED_MODULES configuration option. The value of this configuration option is a comma separated listed of trusted module names. The '*' wildcard can be used at the name of a string to match all strings beginning with the substring before the '*' character. For example 'every*' will make 'every.thing' or 'everything' module trusted. '*' can also be used to make all modules to be trusted. The ".*" wildcard can also be used to match exact modules or submodules names. For example 'every.*' will make 'every' and 'every.thing' modules trusted, but not 'everything'.

Python解释器的链接机制

Currently only CPython 2 and 3 is supported. The GDAL shared object is not explicitly linked at build time to any of the CPython library. When GDAL will need to run Python code, it will first determine if the Python interpreter is loaded in the current process (which is the case if the program is a Python interpreter itself, or if another program, e.g. QGIS, has already loaded the CPython library). Otherwise it will look if the PYTHONSO configuration option is defined. This option can be set to point to the name of the Python library to use, either as a shortname like "libpython2.7.so" if it is accessible through the Linux dynamic loader (so typically in one of the paths in /etc/ld.so.conf or LD_LIBRARY_PATH) or as a full path name like "/usr/lib/x86_64-linux-gnu/libpython2.7.so". The same holds on Windows will shortnames like "python27.dll" if accessible through the PATH or full path names like "c:\python27\python27.dll". If the PYTHONSO configuration option is not defined, it will look for a "python" binary in the directories of the PATH and will try to determine the related shared object (it will retry with "python3" if no "python" has been found). If the above was not successful, then a predefined list of shared objects names will be tried. At the time of writing, the order of versions searched is 2.7, 3.5, 3.6, 3.7, 3.8, 3.9, 3.4, 3.3, 3.2. Enabling debug information (CPL_DEBUG=ON) will show which Python version is used.

即时编译

使用实时编译器可能会显著加快执行时间。 Numba 已成功测试。为了获得更好的性能,建议使用脱机像素函数,以便实时编译器可以缓存其编译。

Given the following mandelbrot.py file :

# Trick for compatibility with and without numba
try:
    from numba import jit
    #print('Using numba')
    g_max_iterations = 100
except:
    class jit(object):
        def __init__(self, nopython = True, nogil = True):
            pass

        def __call__(self, f):
            return f

    #print('Using non-JIT version')
    g_max_iterations = 25

# Use a wrapper for the entry point regarding GDAL, since GDAL cannot access
# the jit decorated function with the expected signature.
def mandelbrot(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                        raster_ysize, r, gt, **kwargs):
    mandelbrot_jit(out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize,
g_max_iterations)

# Will make sure that the code is compiled to pure native code without Python
# fallback.
@jit(nopython=True, nogil=True, cache=True)
def mandelbrot_jit(out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                        raster_ysize, max_iterations):
    ovr_factor_y = float(out_ar.shape[0]) / ysize
    ovr_factor_x = float(out_ar.shape[1]) / xsize
    for j in range( out_ar.shape[0]):
        y0 = 2.0 * (yoff + j / ovr_factor_y) / raster_ysize - 1
        for i in range(out_ar.shape[1]):
            x0 = 3.5 * (xoff + i / ovr_factor_x) / raster_xsize - 2.5
            x = 0.0
            y = 0.0
            x2 = 0.0
            y2 = 0.0
            iteration = 0
            while x2 + y2 < 4 and iteration < max_iterations:
                y = 2*x*y + y0
                x = x2 - y2 + x0
                x2 = x * x
                y2 = y * y
                iteration += 1

            out_ar[j][i] = iteration * 255 / max_iterations

可以使用以下VRT文件(例如用QGIS打开)

<VRTDataset rasterXSize="100000000" rasterYSize="100000000">
    <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
        <PixelFunctionLanguage>Python</PixelFunctionLanguage>
        <PixelFunctionType>mandelbrot.mandelbrot</PixelFunctionType>
        <Metadata>
        <MDI key="STATISTICS_MAXIMUM">255</MDI>
        <MDI key="STATISTICS_MEAN">127</MDI>
        <MDI key="STATISTICS_MINIMUM">0</MDI>
        <MDI key="STATISTICS_STDDEV">127</MDI>
        </Metadata>
        <ColorInterp>Gray</ColorInterp>
        <Histograms>
        <HistItem>
            <HistMin>-0.5</HistMin>
            <HistMax>255.5</HistMax>
            <BucketCount>256</BucketCount>
            <IncludeOutOfRange>0</IncludeOutOfRange>
            <Approximate>1</Approximate>
            <HistCounts>0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|
    0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|
    0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|
    0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|
    0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|
    0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|
    0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0</HistCounts>
        </HistItem>
        </Histograms>
    </VRTRasterBand>
</VRTDataset>

扭曲VRT

扭曲VRT是一个VRTDataset,其子类为“VRTWarpedDataset”。它有一个gdalwarppoptions元素来描述扭曲选项。

<VRTDataset rasterXSize="20" rasterYSize="20" subClass="VRTWarpedDataset">
    <SRS>PROJCS["NAD27 / UTM zone 11N",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26711"]]</SRS>
    <GeoTransform>  4.4072000000000000e+05,  6.0000000000000000e+01,  0.0000000000000000e+00,  3.7513200000000000e+06,  0.0000000000000000e+00, -6.0000000000000000e+01</GeoTransform>
    <Metadata>
        <MDI key="AREA_OR_POINT">Area</MDI>
    </Metadata>
    <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
        <ColorInterp>Gray</ColorInterp>
    </VRTRasterBand>
    <BlockXSize>20</BlockXSize>
    <BlockYSize>20</BlockYSize>
    <GDALWarpOptions>
        <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
        <ResampleAlg>NearestNeighbour</ResampleAlg>
        <WorkingDataType>Byte</WorkingDataType>
        <Option name="INIT_DEST">0</Option>
        <SourceDataset relativeToVRT="1">byte.vrt</SourceDataset>
        <Transformer>
        <ApproxTransformer>
            <MaxError>0.125</MaxError>
            <BaseTransformer>
            <GenImgProjTransformer>
                <SrcGeoTransform>440720,60,0,3751320,0,-60</SrcGeoTransform>
                <SrcInvGeoTransform>-7345.33333333333303,0.0166666666666666664,0,62522,0,-0.0166666666666666664</SrcInvGeoTransform>
                <DstGeoTransform>440720,60,0,3751320,0,-60</DstGeoTransform>
                <DstInvGeoTransform>-7345.33333333333303,0.0166666666666666664,0,62522,0,-0.0166666666666666664</DstInvGeoTransform>
            </GenImgProjTransformer>
            </BaseTransformer>
        </ApproxTransformer>
        </Transformer>
        <BandList>
        <BandMapping src="1" dst="1" />
        </BandList>
    </GDALWarpOptions>
</VRTDataset>

平底锅

2.1 新版功能.

VRT可以描述由 pansharpening operation 泛影VRT将一个全色波段与几个低分辨率的光谱波段结合起来,生成与全色波段相同分辨率的输出光谱波段。

VRT泛影假设全色带和光谱带具有相同的投影(或没有投影)。如果不是这样,则必须在之前的步骤中重新投影。波段可能具有不同的地理转换矩阵,在这种情况下,默认情况下,生成的数据集将具有所有区段的并集作为区段。

目前唯一支持的泛化算法是“加权”Brovey算法。该算法的一般原理是,将谱带重采样到全色带的分辨率后,根据谱带的加权平均值计算出伪全色强度。然后,光谱带的输出值是其输入值乘以实际全色强度与伪全色强度之比。

对应伪码:

pseudo_panchro[pixel] = sum(weight[i] * spectral[pixel][i] for i=0 to nb_spectral_bands-1)
ratio = panchro[pixel] / pseudo_panchro[pixel]
for i=0 to nb_spectral_bands-1:
    output_value[pixel][i] = input_value[pixel][i] * ratio

有效的pansharpen VRT必须将subClass=“vrtpansharpenedataset”声明为VRTDataset top元素的属性。VRTDataset元素必须有一个子元素 PansharpeningOptions 元素。这个PansharpeningOptions元素必须有 PanchroBand 子元素和其中一个 SpectralBand 元素。PanchroBand和spectraband元素必须至少具有 SourceFilename 子元素以指定数据集的名称。他们也可能有 SourceBand 子元素以指定数据集中的带区数(从1开始)。如果未指定,则假定为第一个频带。

spectraband元素通常必须具有 D波段 属性指定输入光谱带必须映射到的输出带的数目(从1开始)。如果未指定该属性,则在计算泛化时将考虑光谱带,但不作为输出带公开。

全色波段和光谱波段通常来自不同的数据集,因为GDAL数据集的波段被假定具有所有相同的维度。光谱带本身可以来自一个或多个数据集。唯一的限制是它们具有相同的尺寸。

下面是一个极简主义工作VRT的例子。它将生成一个数据集,其中3个输出波段对应于multispectral.tif的3个输入光谱波段,并使用pansharped with panchromatic.tif。

<VRTDataset subClass="VRTPansharpenedDataset">
    <PansharpeningOptions>
        <PanchroBand>
            <SourceFilename relativeToVRT="1">panchromatic.tif</SourceFilename>
            <SourceBand>1</SourceBand>
        </PanchroBand>
        <SpectralBand dstBand="1">
            <SourceFilename relativeToVRT="1">multispectral.tif</SourceFilename>
            <SourceBand>1</SourceBand>
        </SpectralBand>
        <SpectralBand dstBand="2">
            <SourceFilename relativeToVRT="1">multispectral.tif</SourceFilename>
            <SourceBand>2</SourceBand>
        </SpectralBand>
        <SpectralBand dstBand="3">
            <SourceFilename relativeToVRT="1">multispectral.tif</SourceFilename>
            <SourceBand>3</SourceBand>
        </SpectralBand>
    </PansharpeningOptions>
</VRTDataset>

在上面的示例中,将从3个声明的输入光谱带创建3个输出pansharpend带。重量是1/3。将使用立方重采样。全色波段的投影和地理变换将被重新用于VRT数据集。

可以创建更明确和声明性的pansharped VRT,例如只允许输出部分输入光谱带(例如,当输入多光谱数据集为RGBNir时,仅允许RGB)。也可以添加“经典”的VRTRasterBands,除了pansharped乐队。

除了上述必需的PanchroBand和SpectralBand元素外,PansharpeningOptions元素还可以具有以下子元素:

  • 算法 :指定平移锐化算法。目前,只支持WeightedBrovey。

  • AlgorithmOptions :指定平移锐化算法的选项。对于WeightedBrovey算法,唯一支持的选项是 砝码 其内容必须是逗号分隔的实值列表的子元素,实值分配每个声明的输入光谱带的权重。必须有与声明的输入光谱带一样多的值。

  • 重采样 :用于将光谱带重采样到全色带分辨率的重采样内核。可以是Cubic(默认)、Average、Near、CubicSpline、Bilinear、Lanczos中的一个。

  • NumThreads: Number of worker threads. Integer number or ALL_CPUS. If this option is not set, the GDAL_NUM_THREADS configuration option will be queried (its value can also be set to an integer or ALL_CPUS)

  • BitDepth :可用于指定全色和光谱带的位深度(例如12)。如果未指定,则将使用全色带中的NBITS元数据项(如果存在)。

  • NoData :全色和光谱波段要考虑的Nodata值。它还将用作输出nodata值。如果未指定,并且所有输入频带都具有相同的nodata值,则将隐式使用该值(除非在nodata中放入特殊的None值以防止出现这种情况)。

  • SpatialExtentAdjustment :可以是 联合 (默认) 交叉NoneNoneWithoutWarning . 控制全色波段和光谱波段具有不同地理空间范围时的行为。默认情况下,并集将采用所有空间范围的并集。交集所有空间范围的交集。“无”将完全不进行任何调整(如果geotransform以某种方式是虚拟的,并且所有标注栏的左上角和右下角都匹配,则可能很有用),但会发出警告。没有警告就等于没有,只是以一种沉默的方式。

下面的示例创建了一个包含4个频带的VRT数据集。第一个波段是全色波段。下面的3个波段是根据具有红、绿、蓝和近红外波段的多光谱栅格计算得到的红、绿、蓝三色盘状波段。计算伪全色强度时考虑了近红外波段,但不受输出波段的限制。

<VRTDataset rasterXSize="800" rasterYSize="400" subClass="VRTPansharpenedDataset">
    <SRS>WGS84</SRS>
    <GeoTransform>-180, 0.45, 0, 90, 0, -0.45</GeoTransform>
    <Metadata>
        <MDI key="DESCRIPTION">Panchromatic band + pan-sharpened red, green and blue bands</MDI>
    </Metadata>
    <VRTRasterBand dataType="Byte" band="1" >
        <SimpleSource>
            <SourceFilename relativeToVRT="1">world_pan.tif</SourceFilename>
            <SourceBand>1</SourceBand>
        </SimpleSource>
    </VRTRasterBand>
    <VRTRasterBand dataType="Byte" band="2" subClass="VRTPansharpenedRasterBand">
        <ColorInterp>Red</ColorInterp>
    </VRTRasterBand>
    <VRTRasterBand dataType="Byte" band="3" subClass="VRTPansharpenedRasterBand">
        <ColorInterp>Green</ColorInterp>
    </VRTRasterBand>
    <VRTRasterBand dataType="Byte" band="4" subClass="VRTPansharpenedRasterBand">
        <ColorInterp>Blue</ColorInterp>
    </VRTRasterBand>
    <BlockXSize>256</BlockXSize>
    <BlockYSize>256</BlockYSize>
    <PansharpeningOptions>
        <Algorithm>WeightedBrovey</Algorithm>
        <AlgorithmOptions>
            <Weights>0.25,0.25,0.25,0.25</Weights>
        </AlgorithmOptions>
        <Resampling>Cubic</Resampling>
        <NumThreads>ALL_CPUS</NumThreads>
        <BitDepth>8</BitDepth>
        <NoData>0</NoData>
        <SpatialExtentAdjustment>Union</SpatialExtentAdjustment>
        <PanchroBand>
            <SourceFilename relativeToVRT="1">world_pan.tif</SourceFilename>
            <SourceBand>1</SourceBand>
        </PanchroBand>
        <SpectralBand dstBand="2">
            <SourceFilename relativeToVRT="1">world_rgbnir.tif</SourceFilename>
            <SourceBand>1</SourceBand>
        </SpectralBand>
        <SpectralBand dstBand="3">
            <SourceFilename relativeToVRT="1">world_rgbnir.tif</SourceFilename>
            <SourceBand>2</SourceBand>
        </SpectralBand>
            <SpectralBand dstBand="4">
            <SourceFilename relativeToVRT="1">world_rgbnir.tif</SourceFilename>
            <SourceBand>3</SourceBand>
        </SpectralBand>
        <SpectralBand> <!-- note the absence of the dstBand attribute, to indicate
                            that the NIR band is not bound to any output band -->
            <SourceFilename relativeToVRT="1">world_rgbnir.tif</SourceFilename>
            <SourceBand>4</SourceBand>
        </SpectralBand>
    </PansharpeningOptions>
</VRTDataset>

多维VRT

3.1 新版功能.

请参阅 多维VRT 页。

vrt:/连接字符串

3.1 新版功能.

在某些上下文中,不必创建文件或提供相当详细的VRT XML内容作为连接字符串,就可以从VRT的特性中获益。为此,自GDAL 3.1以来,数据集名称支持以下URI语法

vrt://{path_to_gdal_dataset}?[bands=num1,...,numN]

例如:

vrt://my.tif?bands=3,2,1

当前唯一受支持的选项是波段。其他可能会在未来增加。

这个选项的作用是改变波段组成。指定的值是源带号(介于1和N之间),可能顺序不正确或重复。这个 mask 值可用于指定全局掩码带区。这也可以看作是运行 gdal_translate -of VRT -b num1 ... -b numN .

多线程问题

警告

下节适用于GDAL<=2.2。从GDAL 2.3开始,VRT数据集的使用遵循标准的GDAL数据集多线程规则(即VRT数据集句柄一次只能由同一线程使用,但可以在同一VRT文件上打开多个数据集句柄并在不同线程中使用它们)

When using VRT datasets in a multi-threading environment, you should be careful to open the VRT dataset by the thread that will use it afterwards. The reason for that is that the VRT dataset uses GDALOpenShared() when opening the underlying datasets. So, if you open twice the same VRT dataset by the same thread, both VRT datasets will share the same handles to the underlying datasets.

SourceFilename上的shared属性指示数据集是否应共享(值为1)或不共享(值为0)。默认值为1。如果在多线程上下文中使用多个引用相同底层源的VRT数据集,则shared应设置为0。或者,可以将VRTU SHAREDU SOURCE configuration(VRTU共享源配置)选项设置为0以强制非共享模式。

性能注意事项

A VRT can reference many (hundreds, thousands, or more) datasets. Due to operating system limitations, and for performance at opening time, it is not reasonable/possible to open them all at the same time. GDAL has a "pool" of datasets opened by VRT files whose maximum limit is 100 by default. When it needs to access a dataset referenced by a VRT, it checks if it is already in the pool of open datasets. If not, when the pool has reached its limit, it closes the least recently used dataset to be able to open the new one. This maximum limit of the pool can be increased by setting the GDAL_MAX_DATASET_POOL_SIZE configuration option to a bigger value. Note that a typical user process on Linux is limited to 1024 simultaneously opened files, and you should let some margin for shared libraries, etc... gdal_translate and gdalwarp, by default, increase the pool size to 450.

驱动程序功能

Supports CreateCopy()

This driver supports the GDALDriver::CreateCopy() operation

Supports Create()

This driver supports the GDALDriver::Create() operation

Supports Georeferencing

This driver supports georeferencing

Supports VirtualIO

This driver supports virtual I/O operations (/vsimem/, etc.)