[GIS] Use gdalbuildvrt then gdal_translate changes Corner Coordinates

coordinatesgdal-translategdalbuildvrtgdalwarp

when I use:

gdalwarp -dstalpha -multi -t_srs epsg:900913 b172.tif b282.tif b770.tif b771.tif b772.tif b832.tif w_25k.tif

I can get one merged tif then go through gdal2tiles for tiling

But if I use gdalwarp it takes really long time to process and very slow.

So I use gdalbuildvrt to create vrt first, then use gdaltranslate to get merged tif.

gdalbuildvrt -overwrite -a_srs epsg:900913 b172.tif b282.tif b770.tif b771.tif b772.tif b832.tif test.vrt

gdal_translate -of GTiff -a_srs epsg:900913 test.vrt  25k.tif

So I get test.tif then go through gdal2tiles the result tif is on top of Africa not the the destination of my tif (Australia).

I use:

gdalinfo 25k.tif

Driver: GTiff/GeoTIFF
Files: 25k.tif
Size is 27063, 13066
Coordinate System is:
PROJCS["Google Maps Global 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.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    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","900913"]]
Origin = (2511424.026000000100000,2339504.385999999900000)
Pixel Size = (0.846642320907504,-0.846706498927300)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 2511424.026, 2339504.386) ( 22d33'37.82"E, 20d33'36.44"N)
Lower Left  ( 2511424.026, 2328441.319) ( 22d33'37.82"E, 20d28' 1.35"N)
Upper Right ( 2534336.707, 2339504.386) ( 22d45'58.80"E, 20d33'36.44"N)
Lower Right ( 2534336.707, 2328441.319) ( 22d45'58.80"E, 20d28' 1.35"N)
Center      ( 2522880.367, 2333972.852) ( 22d39'48.31"E, 20d30'48.92"N)
Band 1 Block=27063x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=27063x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=27063x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=27063x1 Type=Byte, ColorInterp=Alpha

gdalinfo w_25k.tif

Driver: GTiff/GeoTIFF
Files: w_25k.tif
Size is 27080, 13104
Coordinate System is:
PROJCS["Google Maps Global 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.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    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","900913"]]
Origin = (16155902.033114368000000,-4642590.182159157500000)
Pixel Size = (1.080346179541381,-1.080346179541381)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (16155902.033,-4642590.182) (145d 7'51.37"E, 38d26'44.75"S)
Lower Left  (16155902.033,-4656747.038) (145d 7'51.37"E, 38d32'43.07"S)
Upper Right (16185157.808,-4642590.182) (145d23'37.49"E, 38d26'44.75"S)
Lower Right (16185157.808,-4656747.038) (145d23'37.49"E, 38d32'43.07"S)
Center      (16170529.920,-4649668.610) (145d15'44.43"E, 38d29'43.97"S)
Band 1 Block=27080x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=27080x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=27080x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=27080x1 Type=Byte, ColorInterp=Alpha

I found the Corner Coordinates are changes can some one help me to keep the origin Corner Coordinates, please?

enter image description here

The origin Data has *.ers file

DatasetHeader Begin
Version     = "7.1"
Name        = "m734wh.ers"
SourceDataset   = "m734wh.tif"
LastUpdated = Thu May 15 16:07:08 GMT 2014
DataFile    = "m734wh.tif"
DataSetType = Translated
DataType    = Raster
ByteOrder   = LSBFirst
CoordinateSpace Begin
    Datum       = "GDA94"
    Projection  = "LOCAL"
    CoordinateType  = EN
    Rotation    = 0:0:0.0
CoordinateSpace End
RasterInfo Begin
    CellType    = Unsigned8BitInteger
    CellInfo Begin
        Xdimension  = 0.8465509849362
        Ydimension  = 0.846666666666667
    CellInfo End
    NrOfLines   = 1200
    NrOfCellsPerLine    = 1726
    RegistrationCoord Begin
        Eastings    = 2519823.203
        Northings   = 2336158.68
    RegistrationCoord End
    NrOfBands   = 3
    BandId Begin
        Value       = "Red Layer"
    BandId End
    BandId Begin
        Value       = "Green Layer"
    BandId End
    BandId Begin
        Value       = "Blue Layer"
    BandId End
    RegionInfo Begin
        Type        = Polygon
        RegionName  = "All"
        SubRegion   = {
            0   0
            0   1200
            1726    1200
            1726    0
            }
    RegionInfo End
RasterInfo End

DatasetHeader End

Best Answer

Your source data is in GDA94, supposed to be Vicgrid94 covering Philip Island. gdalwarp transforms correctly to a target CRS of epsg:3857, but gdal_translate will not do that.

-a_srs has a different purpose: just assign a CRS, no reprojection.

So you have to use:

gdalwarp -of GTiff -s_srs epsg:3111 -a_srs epsg:3857 test.vrt  25k.tif

BTW I used this batch to merge and tile a set of maps in a local projection:

for %%N in (D:\Karten\gdal\gdal2tiles\NL25\*.tif) DO gdal_translate -of vrt -expand rgba %%N D:\Karten\gdal\gdal2tiles\NL25\%%~nN.vrt
gdalbuildvrt -allow_projection_difference index25.vrt NL25\*.vrt
gdal2tiles --s_srs EPSG:28992 --zoom 15-16 index25.vrt tiles
pause

No big tif file needed, and gdal2tiles reprojects to EPSG:3857 itself. The -expand rgba is needed for paletted source files, so I wrapped those in a vrt for each source file before merging. That might fix your issue in GDAL merge warp has black block can't remove as well, since it adds an alpha channel too.

Related Question