gdal - GDAL2Tiles: MapTiles from BSB/KAP are Switched

As an example I have a BSB file from NOAA which is a navigational map. The map has a SRS of WGS84 and a Mercator projection.

I initially translate the file from BSB format to GTiff format using the following command:

gdal_translate -of GTiff 18554_1.kap noaa.tif

I can view this GTiff file using Paint Shop Pro quite fine. I then translate the file from GTiff format to VRT format using the following command:

gdal_translate -of vrt –expand rgba noaa.tif noaa.vrt

I then try to generate a tile matrix set using the following command:

gdal2tiles.py noaa.vrt tiles_folder

Where the map tiles are located in tiles_folder. If I try to open the openmaps.html or google.html, the map zoom in and out looks good. However, on tiles generator like MBTiles, I see tiles swap. See the image below. I am pulling my hair off to get the tiles to line up. I am guessing it is projection error. I am not understanding those projection types. It would great someone can explain it while finding the error. I am planning to use on a tiling airplane maps/charts app like this. By the way, on GDAL how do I specify the BSB (reference) file that is associated with KAP (data) file? I have seen cases where only one BSB file referred by multiple KAP files.

Answer

I noticed that gdal2tiles numbers the tiles from south to north (according to the TMS specification), while Openstreetmap and others do it from north to south. For my personal use, I changed the code of gdal2tiles to get it right again.

See also: http://osgeo-org.1560.x6.nabble.com/gdal2tiles-tiles-in-wrong-hemisphere-and-or-Openlayers-problem-td3742809.html

---

EDIT

This is what I changed in Version 19288 2010-04-02:

line 1186, insert one line:

ty1=(2**tz - 1) - ty

before

tilefilename = os.path.join(self.output, str(tz), str(tx), "%s.%s" % (ty1, self.tileext))

line 1330, insert

ty2=(2**tz - 1) - ty

before

tilefilename = os.path.join( self.output, str(tz), str(tx), "%s.%s" % (ty2, self.tileext) )

line 1363, insert

y2=(2**(tz+1) - 1) - y

before

dsquerytile = gdal.Open( os.path.join( self.output, str(tz+1), str(x), "%s.%s" % (y2, self.tileext)), gdal.GA_ReadOnly)

Line numbers may have changed in current gdal 1.10.0.