[GIS] Raster Image Tiles created using GDAL2Tiles (from GDAL 2.4.2) does not align properly with Open Street Map Tile layer

gdalgdal2tilesleafletopenstreetmap

I am trying show a custom tile overlay (created from a raster image) on top of openstreetmap base layer using Leaflet.

I used one of the raster images from https://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-1/
and used GDAL2Tiles (I have GDAL 2.4.2.) to create a tile layer. 

gdal2tiles.py -p raster -z 0-8 NE1_HR_LC.tif Natural

The generated tiles are not aligned properly with openstreetmap as base layer. I have attached an image on how it appears when both of the layers appear.Custom tile layer and openstreetmap tile layer

Below is the code in the front end. 

let minzooms =0;
let maxzooms =8;
let beginzooms =3;

var osm = L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {attribution: '&copy; <a href="http://osm.org/copyright">OpenStreetMap</a> contributors', minZoom: minzooms, maxZoom: maxzooms});

var lyr = L.tileLayer('/tiles/{z}/{x}/{-y}.png', {tms: true, opacity: 0.7, attribution: "tiled map", minZoom: minzooms, maxZoom: maxzooms, noWrap: true });

var map = L.map('map', {
        center: [40.86667, 34.56667],
        zoom: beginzooms,
        minZoom: minzooms,
        maxZoom: maxzooms,
        layers: [osm]
    });

 var basemaps = {"OpenStreetMap": osm};
 var overlaymaps = {"Layer": lyr};

 L.control.layers(basemaps, overlaymaps, {collapsed: false}).addTo(map);

 map.fitBounds([[-180, -90], [180, 90]]);

Is gdal2tiles causing this issue?

Best Answer

You are trying to match two different map projections.

The default OpenStreetMap tiles use the EPSG:3857 coordinate reference system, AKA Web Mercator Projection. It looks like this:

mercator projection example

On the other hand, the Natural Earth files use the EPSG:4326 coordianate reference system, AKA equirectangular projection. It looks like this:

equirectangular projection example

The issue of map projections, reference ellipsoids and how those combine to form coordinate reference systems is not trivial if you're new to the GIS field. If you're asking yourself "which is the best projection", the answer is "all of them" (each for a different purpose).

In your particular case, you want gdal2tiles to generate epsg:3857 "mercator" tiles from a epsg:4326 raster source, reprojecting the data during the process.

So, tell gdal2tiles that you want mercator tiles with...

-p mercator

...and tell that the source raster is in EPSG:4326 with...

-s epsg:4326

...and everything together should look like:

gdal2tiles.py -z 0-8 -p mercator -s epsg:4326 NE1_LR_LC.tif Natural

Remember to run gdal2tiles.py --help for more info about command-line parameters.

Related Question