¿Es posible crear VRTWarpedDataset con muchas fuentes usando python y gdal?

Question

Puedo crear fácilmente VRTWarpedDataset para un archivo.

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)

¿Es posible crear un WarpedVRT para muchas fuentes de datos de entrada, como un VRT normal con muchas etiquetas SimpleSource?

He intentado combinar 2 WarpedVRT de tal manera:

<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>

Pero esto no ayudó, el formato XML es incorrecto. La primera figura es VRT me gustaría get.Such VRT no muestra la segunda imagen.

• 

El VRT que me gustaría conseguir

• 

El objetivo de hacer tal VRT:

• VRT simple es demasiado lento para la lectura, debido a BLockSize=128*128px. Y no he encontrado la manera de aumentarlo. Tengo alrededor de 500 dtasources de entrada.
• Evitar los artefactos y la superposición de máscaras entre los límites, cuando se utiliza gdal2tiles y mbutil .

Answer 1

gdalbuildvrt construye un VRT a partir de una lista de conjuntos de datos para hacer mosaicos, así que debería ser lo que buscas.

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

De modo que podría ejecutarlo fácilmente en Python utilizando subprocess o os.system y no sólo (véase Equivalente en Python de gdalbuildvrt ).