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
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
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
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.