2

I've included more information and the scripts that were used to generate the output.

Example results can be found here at the links below and are all generated from the same original GeoTIFF. Each link brings you to the same zoom level and lat/lng in the map. The only difference on the input files is that is that on one of them I slice off all data left of the 180, the other one the map is left complete. The html files are set to the same example point and zoom level to show how that without slicing the file, the image quality is very poor. The one thing I notice is that if I include the 180 meridian, the number of tiles generated is far larger (46k vs 8k).

Bad Results: http://halo3.net/alaska_tiles_no_slice/leaflet.html

Good Results: http://halo3.net/alaska_tiles_sliced_at_180/leaflet.html

Here are the commands I use to generate both tile sets. It will automatically download, unzip and process the input file I am using producing the same results shown at the links above.

How might I be able to get the issue resolved?

I'm still quite a novice at GIS.

#!/bin/bash

wget -SN https://aeronav.faa.gov/content/aeronav/Planning/AK_WallPlan_2.zip && unzip -xo AK_WallPlan_2.zip 
mv "Alaska Wall Planning Chart 2.tif" Alaska_Wall_Planning_Chart_2.tif
gdal_translate  -expand rgb  -strict  -co TILED=YES -co COMPRESS=LZW  -co BIGTIFF=YES Alaska_Wall_Planning_Chart_2.tif alaska_expanded_rgb.vrt
gdalwarp -t_srs EPSG:3857 -r lanczos -overwrite -multi -wo NUM_THREADS=ALL_CPUS -wm 1024 --config GDAL_CACHEMAX 1024 -co TILED=YES alaska_expanded_rgb.vrt alaska_warped.vrt
gdalwarp -t_srs EPSG:3857 -te -19979999 6016546 -13438444 11497579 -r lanczos -overwrite -multi -wo NUM_THREADS=ALL_CPUS -wm 1024 --config GDAL_CACHEMAX 1024 -co TILED=YES alaska_expanded_rgb.vrt alaska_warped_sliced_at_180.vrt
rm -rf alaska_tiles_no_slice
rm -rf alaska_tiles_sliced_at_180
gdalinfo alaska_warped.vrt
time gdal2tiles.py -vvv -z 0-9 --processes=16 -r lanczos -w leaflet alaska_warped.vrt alaska_tiles_no_slice
gdalinfo alaska_warped_sliced_at_180.vrt
time gdal2tiles.py -vvv -z 0-9 --processes=16 -r lanczos -w leaflet alaska_warped_sliced_at_180.vrt alaska_tiles_sliced_at_180

ORIGINAL Post:

I am trying to create a tileset from the following GeoTIFF file.

http://aeronav.faa.gov/content/aeronav/Planning/AK_WallPlan_2.zip

This is an aeronautical chart from the FAA. My issue seems to be that this specific image file does cross the 180th Meridian. This is causing some problems when trying to do the conversion. I'm not sure exactly why, but the resulting tiles have very poor image quality.

converting the image as is

However, if I slice off the image at 180, then it will translate just fine.

Here are the commands I'm using:

gdal_translate -co TILED=TRUE -co COMPRESS=LZW -co BIGTIFF=YES -expand rgba alaska.tif alaska_output.tif
rm alaska_output_translated_slice_at_180.tif
gdalwarp -r lanczos -co LONG_CENTER=90 -s_srs 'PROJCRS["Lambert_Conformal_Conic", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["Lambert Conic Conformal (2SP)", METHOD["Lambert Conic Conformal (2SP)", ID["EPSG",9802]], PARAMETER["Latitude of false origin",60.00305, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8821]], PARAMETER["Longitude of false origin",-141, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8822]], PARAMETER["Latitude of 1st standard parallel",45, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8823]], PARAMETER["Latitude of 2nd standard parallel",65, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8824]], PARAMETER["Easting at false origin",0, LENGTHUNIT["metre",1], ID["EPSG",8826]], PARAMETER["Northing at false origin",0, LENGTHUNIT["metre",1], ID["EPSG",8827]]], CS[Cartesian,2], AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1, ID["EPSG",9001]]], AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1, ID["EPSG",9001]]]]' -t_srs EPSG:3857 -te -19497422 6016546 -13438444 11497579  alaska_output.tif alaska_output_translated_slice_at_180.tif
gdal2tiles.py -z 0-9 --processes=16 -r lanczos -w leaflet -v alaska_output_translated_slice_at_180.tif

rm alaska_output_translated.tif
gdalwarp -r lanczos -co LONG_CENTER=90 -s_srs 'PROJCRS["Lambert_Conformal_Conic", BASEGEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]], CONVERSION["Lambert Conic Conformal (2SP)", METHOD["Lambert Conic Conformal (2SP)", ID["EPSG",9802]], PARAMETER["Latitude of false origin",60.00305, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8821]], PARAMETER["Longitude of false origin",-141, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8822]], PARAMETER["Latitude of 1st standard parallel",45, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8823]], PARAMETER["Latitude of 2nd standard parallel",65, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8824]], PARAMETER["Easting at false origin",0, LENGTHUNIT["metre",1], ID["EPSG",8826]], PARAMETER["Northing at false origin",0, LENGTHUNIT["metre",1], ID["EPSG",8827]]], CS[Cartesian,2], AXIS["easting",east, ORDER[1], LENGTHUNIT["metre",1, ID["EPSG",9001]]], AXIS["northing",north, ORDER[2], LENGTHUNIT["metre",1, ID["EPSG",9001]]]]' -t_srs EPSG:3857 alaska_output.tif alaska_output_translated.tif
gdal2tiles.py -z 0-9 --processes=16 -r lanczos -w leaflet -v alaska_output_translated.tif

map tiled after cutting off 180

Improve this question
1
  • Looks like an interesting error. You could also try creating an issue on github. I don't know how active they are but they might be able to help you there as well. Espescially since you documented how to reproduce the issue github.com/OSGeo/gdal/issues
    – JonasV
    Nov 27, 2020 at 10:39

1 Answer 1

Reset to default
2
+50

Alright, after a bit of testing I think I found the issue. The resolution loss appears in the gdalwarp call. I don't exactly know why the raster is getting written in worse resolution, but it definitely has something to do with how it has to be converted to be displayed in a Pseduo-Mercator projection:

gdalwarp to EPSG:3857

As you can see there's a whole lot of nothing between the eastern and western part. That's also the reason why there's so many more tiles if you don't cut it at 180°.

By specifically setting the target resolution as -tr 169 169, which is the original resolution of your image, you will not get the resolution loss you've seen before. Here's the test code I've used:

gdal_translate -of vrt -expand rgba -mask 1 Alaska_Wall_Planning_Chart_2.tif temp.vrt
gdalwarp -t_srs EPSG:3857 -overwrite -multi -wo NUM_THREADS=ALL_CPUS -wm 1024 -tr 169 169 --config GDAL_CACHEMAX 1024 -co TILED=YES temp.vrt rgb.tif
gdal2tiles -vvv -z 0-9 --processes=8 -r lanczos -w leaflet rgb.tif alaska_tiles_no_slice
Improve this answer
13
  • Thanks. I'm going to try this out, but I think maybe it has something to do with auto-calculating the resolution for such a large area?
    – jr.
    Nov 27, 2020 at 15:18
  • Yes, that's very likely
    – JonasV
    Nov 27, 2020 at 15:19
  • For me, this command is resulting in a 30GB file and it dies at about 10%. I have enough storage. Not sure why it is failing for me.
    – jr.
    Nov 27, 2020 at 22:34
  • Does it die at gdal2tiles?
    – JonasV
    Nov 28, 2020 at 7:25
  • No, outputting the tif in the warp
    – jr.
    Nov 28, 2020 at 11:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.