1

Basically, I have a bunch of different TIFF files, each with their own projection (I think). I want to eventually mosaic and tile them all. But before I tile them I have to combine all the TIFFs in a VRT. When I try to combine them all in to one VRT file, I get the warning gdalbuildvrt does not support heterogeneous projection.

When I try to run gdalbuildvrt with the -allow_projection_difference my resulting VRT file just lays one TIFF on top of the other.

Here is an example output from gdalinfo a.tiff:

Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",46.66666666666666],
    PARAMETER["standard_parallel_2",41.33333333333334],
    PARAMETER["latitude_of_origin",42.16666666666666],
    PARAMETER["central_meridian",-81],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-410068.632339396339376,237081.562232148833573)
Pixel Size = (42.335627217447453,-42.335627217447453)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2017:03:06 15:20:15
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -410068.632,  237081.562) ( 86d 8'15.74"W, 44d11'15.04"N)
Lower Left  ( -410068.632, -287033.503) ( 85d45'32.86"W, 39d28'36.44"N)
Upper Right (  348797.486,  237081.562) ( 76d37'42.21"W, 44d13' 9.57"N)
Lower Right (  348797.486, -287033.503) ( 76d57' 2.60"W, 39d30'22.31"N)
Center      (  -30635.573,  -24975.970) ( 81d22'10.52"W, 41d56'27.88"N)

And output from gdalinfo b.tiff:

Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",46.66666666666666],
    PARAMETER["standard_parallel_2",41.33333333333334],
    PARAMETER["latitude_of_origin",42.166666666666],
    PARAMETER["central_meridian",-89],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-341704.579412317893002,234313.650553743646014)
Pixel Size = (42.335899794430624,-42.335773650774627)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2017:04:04 14:09:40
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -341704.579,  234313.651) ( 93d16'51.86"W, 44d11'51.99"N)
Lower Left  ( -341704.579, -238365.262) ( 92d59'39.89"W, 39d56'47.03"N)
Upper Right (  348370.587,  234313.651) ( 84d38' 8.04"W, 44d11'40.66"N)
Lower Right (  348370.587, -238365.262) ( 84d55'40.03"W, 39d56'36.48"N)
Center      (    3333.004,   -2025.806) ( 88d57'34.76"W, 42d 8'54.28"N)

These two TIFFs portray two regions right next to each other. Since I figured the issue was the the extra parameters (the only thing different about them) I used gdal_translate -a_srs "+proj=lcc +datum=NAD83 +units=m +no_defs to convert them both to the same projection. But this seems to have stripped the georeferencing from the TIFFs.

Here is the output from gdalinfo a.tif after the reprojection:

Coordinate System is:
PROJCS["unnamed",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-410068.632339396339376,237081.562232148833573)
Pixel Size = (42.335627217447453,-42.335627217447453)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2017:03:06 15:20:15
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -410068.632,  237081.562) 
Lower Left  ( -410068.632, -287033.503) 
Upper Right (  348797.486,  237081.562) 
Lower Right (  348797.486, -287033.503) 
Center      (  -30635.573,  -24975.970) 

As you can see, there is no longer any coordinate info attached to each image. So as you would expect, when I then combine them in to a VRT and then render the VRT in QGIS, the two images are once again just stacked on top of each other.

One option that does work is to convert to WGS84 using gdalwarp, but there is a lot of details and typing on my TIFFs that looks really bad and aliased when I reproject to WGS84. I would prefer to keep the TIFFs in LCC and avoid any weird warping.

So how should I go about mosaicing these TIFFs (I have ~30 of them totaling 12 GB) together, when they all have different projections while minimizing distortions?

EDIT: So I found another work flow that is REALLY close to perfect, but not quite there.

  1. Convert PCT tiffs to RGB tiffs: pct2rgb.py <input.tif> <output.tif>
  2. Warp all TIFFs to a uniform projection: `gdalwarp -r lanczos -t_srs EPSG:6703
  3. Build a VRT of all the reprojected TIFFs: gdalbuildvrt <output.vrt> <input1.tif> <input2.tif> ...
  4. Generate tiles: gdal2tiles.py -z <min_zoom>-<max_zoom> -r lanczos <in.vrt> <output_tile_dir>

This produces tiles that are VERY close to my source image.

The source image, zoomed in to show details:The source image, zoomed in to show details

The reprojected TIFF, zoomed in. Note the rotation causes aliasing (most visible along straight lines like the run way indicator in the center):reprojected tiff zoomed in

Unfortunately the compression used by Imgur makes the distortions less visible, but trust me, they are there and pretty glaring on the fully tiled map.

When the map is tiled though, the image is rotated back to its original orientation, but it still has the slight aliasing from the original rotation.

I believe that the cause of these slight distortions is because the image is rotated to adhere to the projection I give it, which I honestly just sort of picked at random. My dataset covers all US states except for Hawaii, and I can't find an LCC projection that covers that wide an area. Is that simply because LCC is not as accurate over a wide area?

Is there some way to combine all these LCC projected TIFFs with different standard parallels, latitudes of origin, and central meridians?

If there is another projection I can convert them to that won't warp the image like WGS84 does (converting to WGS84 smooshed the image vertically), should I be converting to that instead?

Sorry if these are all obvious, and I am just not properly using the tools. I am new to this whole GIS thing =P

  • At some point the VRT will project the images to a single spatial reference, I would think it would be better to take charge of the projection (warping) so it's done in a satisfactory manner and not the arbitrary, default, method. It's not that difficult to warp with a batch file; 12 GB isn't that much, especially if some are already in the correct spatial reference. If you prefer LCC then warp the ones that are not into LCC. Would you like some help writing a batch file or python script to warp? – Michael Stimson Jun 19 '17 at 20:47
  • 1
    You do not re-project with gdal_translate -a_srs "+proj=lcc +datum=NAD83 +units=m +no_defs. You are just assigning wrong srs info into your image. Gdalwarp is the tool for re-projecting. – user30184 Jun 20 '17 at 7:43
  • @MichaelMiles-Stimson they are all in LCC, just with different params (standard parallel, meridian, etc...). The issue is that when I try to warp the images to a single common projection for the VFR, I lose image fidelity. I know this doesn't have to be the case because vfrmap.com somehow managesto combine them all and they look just like the TIFFs given out by the FAA. Check out my edit for more info. – dirtshell Jun 20 '17 at 17:17
  • How about USA Lamberts spatialreference.org/ref/esri/102004 or spatialreference.org/ref/sr-org/7483 to use data in the same spatial reference as Google? Your aliasing comes from the dodgy resample (probably nearest), lossy compression wouldn't help either, use bilinear or cubic to get a better result. Having all your rasters in the same spatial reference as the VRT will speed up refresh if you're serving/consuming in the same spatial reference as your VRT but won't help a bit if you're projecting on the fly to UTM for example. – Michael Stimson Jun 20 '17 at 20:53
  • Changing the resample to cubic did the trick. I was using Lanczos before because I naively assumed it would be best because it was later in the list on the doc. It now looks silky smooth. To be clear, the only reason I am changing the spatial reference is because gdalwarp yells at me about heterogenous projections. But based on my test this should be good. Now I just need to leave my script running over night =) – dirtshell Jun 20 '17 at 22:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.