4

I'm trying to create a seamless mosaic from the FAA's VFR aeronautical charts using gdal_retile.py. The end goal is to serve the mosaic using the ImagePyramid extension in GeoServer. However, I'm seeing gaps as well as a misalignment in the resulting mosaic after running gdal_retile.py, and I'm not sure what I'm doing wrong. I'm still relatively new to GDAL and GIS in general.

To keep things simple, I'm trying to stitch together just 2 rasters (instead of all 54):

Below are the commands I'm using if you want to follow along. I ran these on Ubuntu 14.04 with the gdal-bin and python-gdal packages installed.

1. Download the charts:

This creates a temporary directory structure for all the subsequent commands below, and downloads the original rasters for Cincinnati and Detroit.

mkdir temp && cd temp 
mkdir charts && cd charts 
wget http://aeronav.faa.gov/content/aeronav/sectional_files/Detroit_90.zip 
unzip Detroit_90.zip 
wget http://aeronav.faa.gov/content/aeronav/sectional_files/Cincinnati_94.zip 
unzip Cincinnati_94.zip

2. Download the clipping shapes:

If you look at the original charts, you'll see there are legends around the perimeter. We can't stitch those into the mosaic, we have to clip them out. I have the clipping shapes for the Cincinnati and Detroit charts in my public Dropbox (which is where the commands below download from), but you can alternatively get them from this github repository: aviationCharts (under the clippingShapes/sectional directory; this is where I originally got the clipping shapes from, so credit for those goes to jlmcgraw).

cd .. 
wget https://www.dropbox.com/s/2h2cyzkpih6dxpi/clippingShapes.zip 
unzip clippingShapes.zip

3. Expand colors to RGB:

To create a mosaic, the granules must have a common color model, so this expands both to RGB.

gdal_translate -strict -of VRT -expand rgb "charts/Cincinnati SEC 94.tif" cincinnati_expanded.vrt 
gdal_translate -strict -of VRT -expand rgb "charts/Detroit SEC 90.tif" detroit_expanded.vrt 

4. Crop legends, and reproject to EPSG:3857:

To create a mosaic, the granules must all have a common coordinate system. The original rasters provided by the FAA are each georeferenced to their own local coordinate system (they're all Lambert Conformal Conic, but have different values for the standard parallels and longitude of origin). I will eventually be serving the mosaic using EPSG:3857, so that's what I'm using here.

In addition to reprojecting, this command also uses the clipping shapes downloaded earlier to crop out the legends.

gdalwarp -overwrite -t_srs EPSG:3857 -cutline clippingShapes/Cincinnati_SEC.shp -crop_to_cutline -dstalpha -of vrt -multi -wo NUM_THREADS=ALL_CPUS cincinnati_expanded.vrt cincinnati_warped.vrt 
gdalwarp -overwrite -t_srs EPSG:3857 -cutline clippingShapes/Detroit_SEC.shp -crop_to_cutline -dstalpha -of vrt -multi -wo NUM_THREADS=ALL_CPUS detroit_expanded.vrt detroit_warped.vrt 

5. Create the image pyramid:

This command creates the tiles using the directory/file structure required by the ImagePyramid plugin in GeoServer. This may take a few minutes to run:

mkdir pyramid 
gdal_retile.py -levels 4 -ps 2048 2048 -co "TILED=YES" -co "BLOCKXSIZE=256" -co "BLOCKYSIZE=256" -co "BLOCKYSIZE=256" -s_srs EPSG:3857 -targetDir pyramid cincinnati_warped.vrt detroit_warped.vrt 

Unfortunately, the resulting mosaic has some gaps, and the charts are not lining up correctly. In particular, if you look at the resulting tile pyramid/cincinnati_06_01.tif:

gdal_retile

There's a horizontal gap through the middle (checkered background, which is transparent), and the charts are not lining up (there's an area of overlap between the two charts, and it seems gdal_retile is duplicating instead of just picking one - for example, RANDOLPH CO airport on the left appears twice).

Why is this occurring? gdal_retile is supposed to support overlapping charts, and I think I've satisfied all the other requirements - what am I doing wrong?

If you view cincinnati_warped.vrt and detroit_warped.vrt (the input files I'm feeding to gdal_retile.py in the last step) in QGIS, there's no gap - the images are lining up all right. So, is this an issue with how I'm using gdal_retile, or some sort of limitation that I'm running into? Is there a different command I can use (or adjust my existing command) to get the result I want?

I'm still relatively new to GIS, and especially the GDAL toolkit. The above mentioned github project (aviationCharts) has helped me a lot with figuring this out, but I'm trying to take a slightly different approach (GeoServer + ImagePyramid, instead of mbtiles), and also trying to understand the process and figure out why things aren't working the way I thought.

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.