Virtual Raster format - VRT¶
Mezi mnoha formáty, které podporuje knihovna GDAL má Virtual rastr driver zvláštní místo. Jedná se soubor XML, kterým můžeme zpracovávat vstupní rastrové soubory za běhu. Soubor tak nezabírá místo. VRT soubory lze řetězit, přidávat do nich kód v jazyce Python a tak podobně.
Transformace souřadnicového systému¶
Jako první příklad prostě pouze ztransformujeme souřadnicový systém vstupního souboru:
gdalwarp -t_srs epsg:5514 -of VRT _data/Copernicus_DSM_COG_30_N49_00_E015_00_DEM.tif _data/Copernicus_DSM_COG_30_N49_00_E015_00_DEM-krovak.vrt
gdalinfo _data/Copernicus_DSM_COG_30_N49_00_E015_00_DEM-krovak.vrt
Podíváme-li se na VRT soubor blíže, uvidíme, že se jedná o celkem malý soubor XML
<VRTDataset rasterXSize="1093" rasterYSize="1524" subClass="VRTWarpedDataset">
<SRS dataAxisToSRSAxisMapping="1,2">PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_of_the_Unified_Trigonometrical_Cadastral_Network",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6156"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.8333333333333],PARAMETER["azimuth",30.2881397527778],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5514"]]</SRS>
<GeoTransform> -7.1734568408439704e+05, 7.8201966936414763e+01, 0.0000000000000000e+00, -1.0580974742596082e+06, 0.0000000000000000e+00, -7.8201966936414763e+01</GeoTransform>
<Metadata>
<MDI key="AREA_OR_POINT">Point</MDI>
</Metadata>
<VRTRasterBand dataType="Float32" band="1" subClass="VRTWarpedRasterBand">
<ColorInterp>Gray</ColorInterp>
</VRTRasterBand>
<BlockXSize>512</BlockXSize>
<BlockYSize>128</BlockYSize>
<GDALWarpOptions>
<WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
<ResampleAlg>NearestNeighbour</ResampleAlg>
<WorkingDataType>Float32</WorkingDataType>
<Option name="INIT_DEST">0</Option>
<Option name="ERROR_OUT_IF_EMPTY_SOURCE_WINDOW">FALSE</Option>
<SourceDataset relativeToVRT="1">Copernicus_DSM_COG_30_N49_00_E015_00_DEM.tif</SourceDataset>
<Transformer>
<ApproxTransformer>
<MaxError>0.125</MaxError>
<BaseTransformer>
<GenImgProjTransformer>
<SrcGeoTransform>14.9995833333333337,0.000833333333333333387,0,50.0004166666666663,0,-0.000833333333333333387</SrcGeoTransform>
<SrcInvGeoTransform>-17999.5,1200,0,60000.4999999999927,0,-1200</SrcInvGeoTransform>
<DstGeoTransform>-717345.684084397042,78.2019669364147632,0,-1058097.47425960819,0,-78.2019669364147632</DstGeoTransform>
<DstInvGeoTransform>9172.98774169789976,0.012787402148248906,0,-13530.3179154040554,0,-0.012787402148248906</DstInvGeoTransform>
<ReprojectTransformer>
<ReprojectionTransformer>
<SourceSRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SourceSRS>
<TargetSRS>PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_of_the_Unified_Trigonometrical_Cadastral_Network",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6156"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.8333333333333],PARAMETER["azimuth",30.2881397527778],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5514"]]</TargetSRS>
<Options>
<Option key="CENTER_LONG">15.4996</Option>
<Option key="AREA_OF_INTEREST">14.99958333333333,49.00041666666667,15.99958333333333,50.00041666666667</Option>
</Options>
</ReprojectionTransformer>
</ReprojectTransformer>
</GenImgProjTransformer>
</BaseTransformer>
</ApproxTransformer>
</Transformer>
<BandList>
<BandMapping src="1" dst="1" />
</BandList>
</GDALWarpOptions>
</VRTDataset>
Soubor lze bez problému načíst do QGIS a zobrazit.
Výhodou je malá velikost. Nevýhodou pomalejší zpracování
Spojení více souborů do jednoho - gdalbuildvrt¶
Pomocí gdalbuildvrt můžeme spojit více souborů do jednoho
gdalbuildvrt _data/dem.vrt _data/Copernicus_DSM_COG_30_N49_00_E015_00_DEM.tif _data/Copernicus_DSM_COG_30_N50_00_E015_00_DEM.tif
gdalinfo _data/dem.vrt
Vložení kódu Python do VRT¶
Do VRT souboru lze přímo vložit kód v jazyce Python. Kód lze i udržovat v externím souboru. Příklad je převod výšky v metrech na stopy:
<VRTDataset rasterXSize="1200" rasterYSize="1200">
<SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
<GeoTransform> 1.4999583333333334e+01, 8.3333333333333339e-04, 0.0000000000000000e+00, 5.0000416666666666e+01, 0.0000000000000000e+00, -8.3333333333333339e-04</GeoTransform>
<VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
<PixelFunctionType>feets</PixelFunctionType>
<PixelFunctionLanguage>Python</PixelFunctionLanguage>
<PixelFunctionCode><![CDATA[
import numpy as np
def feets(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
out_ar[:] = in_ar[0]*3.28084
]]>
</PixelFunctionCode>
<SimpleSource>
<SourceFilename relativeToVRT="1">Copernicus_DSM_COG_30_N49_00_E015_00_DEM.tif</SourceFilename>
</SimpleSource>
</VRTRasterBand>
</VRTDataset>
A můžeme se pokusit získat statistické údaje
gdalinfo _data/Copernicus_DSM_COG_30_N49_00_E015_00_DEM-feets.vrt -stats
[...]
Metadata:
STATISTICS_MAXIMUM=2752.158203125
STATISTICS_MEAN=1673.4303891539
STATISTICS_MINIMUM=648.65942382812
STATISTICS_STDDEV=365.32895543763
STATISTICS_VALID_PERCENT=100
Varování
Narazíme na problém, že kvůli bezpečnosti nelze vykonávat kód v Pythonu jen tak. Je potřeba nastavit proměnnou GDAL_VRT_ENABLE_PYTHON na hodnotu YES. Na os Linux: export GDAL_VRT_ENABLE_PYTHON=YES
Soubor by měl jít zobrazit i v QGIS.