1

I have warped some images from various state plane coordinate to EPSG 3857, but the output of gdalbuildvrt is messed up (see image).

I've tried dozens of different combinations of gdalwarped and gdalbuiltvrt.

If I gdalwarp them without -dstalpha and use -srcnodata "255 255 255", then the gdalbuiltvrt works. The problem with this method is that there are some solid white pixels that I don't want as no data.

Any ideas why the VRT is bugged out when dstalpha is used?

gdalwarp -t_srs EPSG:3857 INPUT.tif OUTPUT.tif -dstalpha
gdalbuildvrt -addalpha -overwrite OUTPUT.vrt INPUT.tif INPUT2.tif

The gdalinfo on these 

Files: MERGEDVRT.vrt
       INPUT1.tif
       INPUT2.tif
       INPUT3.tif
       INPUT4.tif
Size is 12447, 14841
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 = (-10715103.831728048622608,4113609.510417979210615)
Pixel Size = (1.215874674612295,-1.215874674612295)
Corner Coordinates:
Upper Left  (-10715103.832, 4113609.510) ( 96d15'19.50"W, 34d37'45.26"N)
Lower Left  (-10715103.832, 4095564.714) ( 96d15'19.50"W, 34d29'44.69"N)
Upper Right (-10699969.840, 4113609.510) ( 96d 7'10.07"W, 34d37'45.26"N)
Lower Right (-10699969.840, 4095564.714) ( 96d 7'10.07"W, 34d29'44.69"N)
Center      (-10707536.836, 4104587.112) ( 96d11'14.78"W, 34d33'45.07"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  Min=0.000 Max=240.000
  Minimum=0.000, Maximum=240.000, Mean=105.003, StdDev=38.946
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=240
    STATISTICS_MEAN=105.00307761725
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=38.945866179465
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  Min=0.000 Max=251.000
  Minimum=0.000, Maximum=251.000, Mean=112.130, StdDev=37.466
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=251
    STATISTICS_MEAN=112.12989263679
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=37.465593515256
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  Min=0.000 Max=239.000
  Minimum=0.000, Maximum=239.000, Mean=85.818, StdDev=29.318
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=239
    STATISTICS_MEAN=85.818457335167
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=29.317767453532
Band 4 Block=128x128 Type=Byte, ColorInterp=Alpha






Files: INPUTX.tif
Size is 6722, 7885
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 = (-10715103.124073004350066,4105152.986530220136046)
Pixel Size = (1.215441657313418,-1.215441657313418)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_ARTIST=Surdex Corporation, 636-368-4400, www.surdex.com
  TIFFTAG_DATETIME=2017:11:07 09:27:51
  TIFFTAG_DOCUMENTNAME=Coalgate SW 3409631
  TIFFTAG_HOSTCOMPUTER=DPNODE8
  TIFFTAG_IMAGEDESCRIPTION=Image courtesy of USDA Farm Service Agency's National Agriculture Imagery Program (NAIP) under Aerial Photography Field Office contract AG-8447-C-12-0003. Imagery has been placed in the public domain and may be used and reproduced without permission or fee. Please credit 'NAIP imagery provided by USDA Farm Service Agency' on any use.
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-10715103.124, 4105152.987) ( 96d15'19.47"W, 34d34' 0.14"N)
Lower Left  (-10715103.124, 4095569.229) ( 96d15'19.47"W, 34d29'44.81"N)
Upper Right (-10706932.925, 4105152.987) ( 96d10'55.25"W, 34d34' 0.14"N)
Lower Right (-10706932.925, 4095569.229) ( 96d10'55.25"W, 34d29'44.81"N)
Center      (-10711018.025, 4100361.108) ( 96d13' 7.36"W, 34d31'52.50"N)
Band 1 Block=6722x1 Type=Byte, ColorInterp=Red
  Min=0.000 Max=252.000
  Minimum=0.000, Maximum=252.000, Mean=105.288, StdDev=47.312
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=252
    STATISTICS_MEAN=105.28834341631
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=47.311851520036
Band 2 Block=6722x1 Type=Byte, ColorInterp=Green
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=110.606, StdDev=45.391
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=110.60634566432
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=45.391333364765
Band 3 Block=6722x1 Type=Byte, ColorInterp=Blue
  Min=0.000 Max=255.000
  Minimum=0.000, Maximum=255.000, Mean=87.709, StdDev=37.475
  Mask Flags: PER_DATASET ALPHA
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=87.709082944891
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=37.474871475178
Band 4 Block=6722x1 Type=Byte, ColorInterp=Alpha

EDIT: "Source Images" meaning the output of gdalwarp. The original images are in SPC, which are listed at the bottom of the post.

enter image description here

ORIGINAL IMAGES:

http://s3.amazonaws.com/pian/original1.tif

http://s3.amazonaws.com/pian/original2.tif

http://s3.amazonaws.com/pian/original3.tif

http://s3.amazonaws.com/pian/original4.tif

Improve this question
3
  • Where you aim to use your images? Do you really need the physically warped EPSG:3857 tiff versions? If you use your images in QGIS you could do first gdalbuildvrt temp.vrt *.tif and then gdalwarp -t_srs EPSG:3857 temp.vrt final.vrt. Commented Apr 5, 2019 at 16:41
  • Yes I'm processing it to be on a web map, and EPSG 3857 is required. Strangely, the EPSG didn't seem to be the issue, but only the difference of -dstalpha and -srcnodata "255 255 255" (this works, but has other problems) Commented Apr 5, 2019 at 16:58
  • I would write mail to gdal-dev mailing list and ask if this could mean a bug. As far as I understand it should work even without -addalpha, and band 4 is correctly marked as Alpha in the VRT file. Materializing the VRT into geotiff with gdal_translate gives similar result. Commented Apr 5, 2019 at 20:31

1 Answer 1

Reset to default
0

I think the issue here is that gdalbuildvrt does not do alpha compositing. Each of your input images have transparent areas around the edges from the reprojection step. In the resulting mosaic, the transparent areas in the top image will override the actual data in the images underneath.

Unless VRT is the only format you can use, you could do the entire thing in one process (example with GeoPackage as output): 

gdalwarp -of GPKG path\to\sourceimages\*.tif output.gpkg

In my experience, gdalwarp properly handles multiple overlapping images with alpha bands.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.