Is it possible to create VRTWarpedDataset with many sources using python and gdal?

Question

I can easely create VRTWarpedDataset for one file.

import gdal,osr
in_ds = gdal.Open('/path/to/input',gdal.GA_ReadOnly)
out_drv = gdal.GetDriverByName('VRT')
prj = in_ds.GetProjection()
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326) 
warped = gdal.AutoCreateWarpedVRT(in_ds,prj,sr.ExportToWKT())
out_ds = out_drv.CreateCopy('/path/to/out.vrt',warped)

Is it possible to create WarpedVRT for many input datasources, like usual VRT with many SimpleSource tags??

i've tried to combine 2 WarpedVRT in such way:

<VRTDataset rasterXSize="8192" rasterYSize="16384" subClass="VRTWarpedDataset">
 <SRS>PROJCS["WGS 84 / UTM zone 37N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",39],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","32637"]]</SRS>
 <GeoTransform> 6.4417049775006599e+005, 4.9999999999975835e-002, 0.0000000000000000e+000, 6.0496980602764217e+006, 0.0000000000000000e+000,-4.9999999999975835e-002</GeoTransform>
 <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
 <ColorInterp>Red</ColorInterp>
 </VRTRasterBand>
 <VRTRasterBand dataType="Byte" band="2" subClass="VRTWarpedRasterBand">
 <ColorInterp>Green</ColorInterp>
 </VRTRasterBand>
 <VRTRasterBand dataType="Byte" band="3" subClass="VRTWarpedRasterBand">
 <ColorInterp>Blue</ColorInterp>
 </VRTRasterBand>
 <BlockXSize>8192</BlockXSize>
 <BlockYSize>128</BlockYSize>
 <GDALWarpOptions>
 <WarpMemoryLimit>6.71089e+007</WarpMemoryLimit>
 <ResampleAlg>NearestNeighbour</ResampleAlg>
 <WorkingDataType>Byte</WorkingDataType>
 <Option name="INIT_DEST">0</Option>
 <SourceDataset relativeToVRT="0">path/to/first.jpg</SourceDataset> 
 <Transformer>
 <GenImgProjTransformer>
 <SrcGeoTransform>644170.49775006599,0.050000000000000003,0,6049698.0602764217,0,-0.050000000000000003</SrcGeoTransform>
 <SrcInvGeoTransform>-12883409.955001319,20,0,120993961.20552842,0,-20</SrcInvGeoTransform>
 <DstGeoTransform>644170.49775006599,0.049999999999975835,0,6049698.0602764217,0,-0.049999999999975835</DstGeoTransform>
 <DstInvGeoTransform>-12883409.955007546,20.000000000009667,0,120993961.20558691,0,-20.000000000009667</DstInvGeoTransform>
 </GenImgProjTransformer>
 </Transformer>
 <BandList>
 <BandMapping src="1" dst="1"/>
 <BandMapping src="2" dst="2"/>
 <BandMapping src="3" dst="3"/>
 </BandList>
 </GDALWarpOptions>
 <GDALWarpOptions>
 <WarpMemoryLimit>6.71089e+007</WarpMemoryLimit>
 <ResampleAlg>NearestNeighbour</ResampleAlg>
 <WorkingDataType>Byte</WorkingDataType>
 <Option name="INIT_DEST">0</Option>
 <SourceDataset relativeToVRT="0">path/to/second.jpg</SourceDataset>
 <Transformer>
 <GenImgProjTransformer>
 <SrcGeoTransform>644170.49775006599,0.050000000000000003,0,6049288.4602764221,0,-0.050000000000000003</SrcGeoTransform>
 <SrcInvGeoTransform>-12883409.955001319,20,0,120985769.20552844,0,-20</SrcInvGeoTransform>
 <DstGeoTransform>644170.49775006599,0.049999999999975835,0,6049288.4602764221,0,-0.049999999999975835</DstGeoTransform>
 <DstInvGeoTransform>-12883409.955007546,20.000000000009667,0,120985769.20558691,0,-20.000000000009667</DstInvGeoTransform>
 </GenImgProjTransformer>
 </Transformer>
 <BandList>
 <BandMapping src="1" dst="1"/>
 <BandMapping src="2" dst="2"/>
 <BandMapping src="3" dst="3"/>
 </BandList>
 </GDALWarpOptions>
</VRTDataset>

But this didn't helped, XML format is wrong. The first figure is VRT i would like to get.Such VRT doesn't show the second image.

• VRT warped with two datasources

The VRT i would like to get

• The VRT i would like to get

The aim for doing such VRT:

• Simple VRT is too slow for reading, because of BLockSize=128*128px. And i didn't find way to increase it. I have about 500 input dtasources.
• Avoid of artifacts and mask overlay between boundaries, when using gdal2tiles and mbutil.

Answer 1

gdalbuildvrt builds a VRT from a list of datasets in order to mosaic them, so it should be what you're looking for.

gdalbuildvrt /path/to/mosaic.vrt /path/to/out*.vrt

So you could easily execute it in Python using subprocess or os.system and not only (see Python equivalent of gdalbuildvrt).

• Dmitriy Litvinov: I've tried to do this, but It will create simple VRT with <SimpleSource>/path/to/warped.vrt</SimpleSource> sources, where warped vrt links to original raster (tiff,jpeg,etc). But the resulting VRT is still of 128*128px BlockSize and is still too slow for reading.Mayby there is another way to read lots of rasters as\like chunk or unified raster?
• Antonio Falciano: It seems that gdalbuildvrt considers the warped VRTs as <SimpleSource> correctly, because there's not rescaling and/or offsetting the range of the source values. Instead, to change the block size you can use gdal_translate with the GeoTIFF creation options (-co) BLOCKXSIZE and BLOCKYSIZE.
• Dmitriy Litvinov: Gdal_translate will create a copy of my data, it will not use rational of disk space. I have many folders with about 500 8192 * 8192px inputs in each folder. Ideally i would like to improove gdal2tiles script to have many files as input. VRT is a good variant, but low perfomance (WarpedVRT has good perfomance but one file as input). When i use gdal2tiles for each file independantly and then merge them the output files have artifacts and overlay mask between boundaries. Maybe there is way to stream reading (or scanaline of neighbour files) of many inputs?