Create Blue Marble XYZ tiles

Question

I am trying to create my own XYZ tiles to use them with OpenLayers 3 as my local tiles. On this example you can see the objective/idea.

I want to use the NASA Blue Marble Earth Visible maps, as an example this one from July 2004

On the readme PDF it says that the image is on projection datum WGS84. So the as we know this projection bounds are on EPSG:4326.

I georeference the image on this way.

gdal_translate -of GTiff -a_srs EPSG:4326 -a_ullr -180 90 180 -90 world.topo.bathy.200407.3x5400x2700.jpg world.topo.bathy.200407.3x5400x2700.tif

The information of the new generated GeoTIFF can be checked with:

$ gdalinfo world.topo.bathy.200407.3x5400x2700.tif
Driver: GTiff/GeoTIFF
Files: world.topo.bathy.200407.3x5400x2700.tif
Size is 5400, 2700
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.066666666666667,-0.066666666666667)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=5400x1 Type=Byte, ColorInterp=Red
Band 2 Block=5400x1 Type=Byte, ColorInterp=Green
Band 3 Block=5400x1 Type=Byte, ColorInterp=Blue

As I want to create XYZ tiles with OpenLayers 3, checking the documentation, I see I need to warp the projection from EPSG:4326 to EPSG:3857, which I do with:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 world.topo.bathy.200407.3x5400x2700.tif Marble_3857.tif

The problem is that the shape of this final image is (4009 x 4514) so it is not square which should be for XYZ tiles.

I have seen from MapTiler this example that uses that map which is already tiled in 256px tiles, something similar to what I want.

First of all, I have no square shape image projection so I cannot make tiles. And secondly, how can I make different zoom level tiles (zoom till 8 level)? Is there any tool on GDAL for that? I tried gdal2tiles but does not work with following command:

gdal2tiles.py --s_srs EPSG:3857 --zoom 0-3 -p 'raster' Marble_3857.tif

EDIT

Following the answer of AndreJ it mostly seems to work. First I used gdal_translate as mentioned:

gdal_translate -of GTiff -a_srs EPSG:4326 -a_ullr -180 90 180 -90 world.topo.bathy.200407.3x5400x2700.jpg world.topo.bathy.200407.3x5400x2700.tif

Then, before gdal2tiles, as I am using Ubuntu I had to change something on the python code as mentioned here to avoid errors on projection.

After that, I could run the command of AndreJ:

gdal2tiles.py --s_srs EPSG:4326 --zoom 0-5 world.topo.bathy.200407.3x5400x2700.tif ./tiles

But the result is quite strange because it leaves a transparency column on the right on all the zoom levels... (can be checked downloading the below picture)

Which can be the reason? The modified code? Any option I should define?

Answer 1

There is no need to reproject the raster. This one worked for me:

gdal2tiles --s_srs EPSG:4326 --zoom 0-5 bluemarble.tif D:\Tiles\bluemrble

gdal2tiles automatically cuts off the pole areas at 85.0511 degrees to get a square output at zoom level 0:

Answer 2

I have prepared ESRI World files for the Nasa Blue marble dataset: https://gist.github.com/klokan/b654b18367f372d58c58

If you place them next to the downloaded files they will get correct geolocation automatically.

MapTiler (http://www.maptiler.com) will generate the correct map tiles for you directly after you drag&drop the files into the user interface. It has also available command line interface for automation and installers for all different platforms like Windows, Mac OS X, various Linux distributions and Docker.

A free demo of MapTiler Pro is available at http://www.maptiler.com/demo/.

MapTiler is most efficient when you supply the individual input files directly - so do not make any preprocessing or merging in advance.

If you want to use open-source GDAL library tools you would probably first need to merge the individual input files into one VRT (with gdalbuildvrt), then reproject correctly (with gdalwarp) into the exact extent of the Spherical Mercator, then tile with my student project gdal2tiles.py. The result would be png or jpg tiles which may need size optimisation with pngnq or similar utility.

MapTiler is doing internally similar steps, but more efficiently and significantly faster. It also automatically handle dateline crossing datasets, is able to run on cluster and is easy to install and use and regularly updated. Try it and you will see! ;-)

Disclaimer: I am the original author of gdal2tiles.py, MapTiler and I have also rendered the tiles used in the jsfiddle example mentioned in the question.