Raster extents incorrect for vrt when changing projection with -a_srs flag

Question

I want to generate a .vrt in epsg:3857 from a series of raster tiles in epsg:4326, and am trying to do so in one step using gdalbuildvrt with the -a_srs flag:

gdalbuildvrt -a_srs epsg:3857 ortho_3857.vrt ortho-*.tif

The projection is changed to 3857 as desired, however the corner coordinates, center and pixel sizes are incorrect (they are still in 4326, with the geographic extents of the corner coords messed up):

Driver: VRT/Virtual Raster
Files: ortho_3857.vrt
       /Users/jeremyeastwood/Desktop/debug/ortho-0-0.tif
       /Users/jeremyeastwood/Desktop/debug/ortho-0-1.tif
       /Users/jeremyeastwood/Desktop/debug/ortho-1-0.tif
       /Users/jeremyeastwood/Desktop/debug/ortho-1-1.tif
Size is 9191, 7284
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    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["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-122.318962274868440,37.873466153295041)
Pixel Size = (0.000000446142119,-0.000000446142119)
Corner Coordinates:
Upper Left  (-122.3189623,  37.8734662) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Lower Left  (-122.3189623,  37.8702165) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Upper Right (-122.3148618,  37.8734662) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Lower Right (-122.3148618,  37.8702165) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Center      (-122.3169120,  37.8718413) (  0d 0' 3.96"W,  0d 0' 1.22"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  NoData Value=0
Band 4 Block=128x128 Type=Byte, ColorInterp=Alpha
  NoData Value=0 

I can achieve the desired result by generating the vrt in the native projection then warping the result to 3857 with gdalwarp afterwards, but then I need to keep both .vrts about for future operations as the warped vrt references the native one, which references the original geotiffs, which is pretty messy.

Any suggestions on correct usage of the -a_srs flag in gdalbuildvrt (or a better understanding of its limitations), or how to generate a single warped vrt in a different projection from the raw geotiffs?

Answer 1

-a_srs just assigns the given SRS into metadata but nothing else.

I believe this is what you did already and I don't know any better way for doing it:

gdalbuildvrt -a_srs epsg:4326 ortho_4326.vrt ortho-*.tif

gdalwarp -of VRT -s_srs epsg:4326 -t_srs epsg:3857 ortho_4326.vrt ortho_3857.vrt.

There is one exception, if your final goal is to use VRT as a mosaic for MapServer you can achieve what you want with tileindex system by using gdaltindex with the -t_srs switch http://www.gdal.org/gdaltindex.html.