Fixed geodetic tile calculation and cleaned up code.

[stevendlander]

The geodetic profile was creating tile caches with a resolution factor of 0.703125, instead of the standard for TMS geodetic tiles which is 1.40625. This error broke tile numbering for subsequent zoom levels.

Additionally, refactored code within the geodetic object for better maintainability.

[rouault]

It is not obvious to me that this patch is needed.

I have just generated the following :
python swig/python/scripts/gdal2tiles.py -p geodetic world_4326.tif gdal2tiles_output

And when I display gdal2tiles_output/openlayers.html my GeoTIFF is correctly overlayed onto the VMAP layer (after replacing labs.metacarta.com by tilecache.osgeo.org. See http://trac.osgeo.org/gdal/changeset/26162 )

[stevendlander]

The openlayers geodetic code handles this on line 2207 by subtracting 1 from the server zoom level. This is not necessary, and removing that would make it mirror the mercator code on line 2183.

I would update the pull request but I don't know how, now that it is submitted. :(
EDIT: My ignorance has been overcome. :)

[stevendlander]

Added a fix for the openlayers browser to the pull request.

A little backstory:

I have been working with MapProxy extensively as a way to turn tiles into sqlite databases. Most of my data is in WGS84 global geodetic format natively, yet in the past I was tiling this data to mercator and then back to WGS84 when creating the sqlite dataset. When I attempted to tile to geodetic directly via gdal2tiles, the result did not match the default profile for global geodetic:

Configuration:
bbox_: [-180.0, -90.0, 180.0, 90.0]
origin_: 'sw'
srs: 'EPSG:4326'
tile_size*: [256, 256]
Levels: Resolutions, # x * y = total tiles
00: 1.40625, # 1 * 1 = 1
01: 0.703125, # 2 * 1 = 2
02: 0.3515625, # 4 * 2 = 8
03: 0.17578125, # 8 * 4 = 32
...
18: 5.364418029785156e-06, # 262144 * 131072 = 34.360G
19: 2.682209014892578e-06, # 524288 * 262144 = 137.439G

Instead, gdal2tiles is calling the level with 2 tiles as level 0. The leaflet group has experienced a similar issue with how MetaCarta handles TMS here: Leaflet/Leaflet#1207 (comment).

You can also see here that OpenLayers make this calculation with 1.40625 being the default value: . While setting the maxResolution field at 0.703125 would indeed fix this issue, I felt it made sense to adhere to a more widely accepted resolution factor.

[rouault]

Committed in r26164 and r26165

[rouault]

Steven,

there's one aspect that has been overlooked. It is related to generating TMS tiles. What was generated previously was compatible (if one excepts the existing bugs in the x/y inversion in the generated tilemapresource.xml ) of a global-geodetic profile ( http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification#global-geodetic ). This is no longer the case.

[stevendlander]

Understood, but to my knowledge TMS is not an official OGC standard and thus could have many interpretations. OSGeo has made "level 0 = 0.703125" their interpretation for compatibility with their products.

I was able to find documentation from OGC regarding the successor to TMS, which is WMTS. That specification calls for the resolution change I have outlined here. See page 119 of the PDF located here: http://www.opengeospatial.org/standards/wmts, titled "OpenGIS Web Map Tile Service Implementation Standard."

So will you be keeping the implementation as-is? From the way this discussion has evolved I think either way could be considered correct.

[rouault]

I'm wondering if there shouldn't be a flag for the user to specifically ask for TMS compatibility, in which case the old behaviour would be kept. Otherwise we would adopt the new one.

[stevendlander]

I will update code to support such a flag.

[rouault]

The 2 above commits do not work :

python /home/even/osgeo-gdal-git-full/gdal/swig/python/scripts/gdal2tiles.py -p geodetic -d world_4326.tif world_4326_tms
Traceback (most recent call last):
File "/home/even/osgeo-gdal-git-full/gdal/swig/python/scripts/gdal2tiles.py", line 2277, in
gdal2tiles.process()
File "/home/even/osgeo-gdal-git-full/gdal/swig/python/scripts/gdal2tiles.py", line 481, in process
self.open_input()
File "/home/even/osgeo-gdal-git-full/gdal/swig/python/scripts/gdal2tiles.py", line 994, in open_input
self.geodetic = GlobalGeodetic() # from globalmaptiles.py
File "/home/even/osgeo-gdal-git-full/gdal/swig/python/scripts/gdal2tiles.py", line 358, in init
if self.options.tmscompatible:
AttributeError: 'GlobalGeodetic' object has no attribute 'options'

[stevendlander]

Yes I am aware, I am not done yet. I was hoping to do it from the git editor but it will take a bit more work on my end. Should be finished up soon.

[stevendlander]

Should be good to go now.

[rouault]

Merged 1207814, 9f37a9d and 36af00c in r26204 . I'm confused about b007343 and 2a6b95f ... Are we OK or are there still commits to come ?

[stevendlander]

All done, thanks for working with me on this. :)

[rouault]

The code in GlobalGeodetic() is a bit confusing. In the "if tmscompatible is not None:" branch, we use non-TMS resFact, whereas in "if tmscompatible is None" branch we use TMS settings.

[stevendlander]

I coded the flag so that if the user adds -d, they will get the behavior that the script currently produces. If -d is omitted, then these changes will go in effect.

Would you like that order swapped?

[rouault]

I'm just saying that the code is confusing. I think the comments are inverted in the GlobalGeodetic() constructor.

That should likely be :
if tmscompatible is not None:
# Defaults the resolution factor to 0.703125 (2 tiles @ level 0)
# Adhers to OSGeo TMS spec http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification#global-geodetic
self.resFact = 180.0 / self.tileSize
else:
# Defaults the resolution factor to 1.40625 (1 tile @ level 0)
# Adheres OpenLayers, MapProxy, etc default resolution for WMTS
self.resFact = 360.0 / self.tileSize

[stevendlander]

Yes, the comments were backwards. I have been awake far too long. :)

[rouault]

ok, commited in r26205