Search code examples
d3.jsgeojsonrastergdal

Aligning Natural Earth Geojson and Raster to render in D3

I am trying to render the world map with elevation data using D3.

For this I use Natural Earth 50m land geojson : https://github.com/martynafford/natural-earth-geojson/tree/master/50m/physical

And Natural Earth elevation raster data : https://www.naturalearthdata.com/downloads/50m-raster-data/50m-shaded-relief/

I am using this tutorial : https://datawanderings.com/2020/08/08/raster-backgrounds/

So I first found the bounds of the geojson :

ogrinfo ne_50m_land.json -so -al | grep Extent  
Extent: (-180.000000, -89.998926) - (180.000000, 83.599609)

I then cut the TIF file :

gdal_translate -projwin -180.000000 83.599609 180.000000 -89.998926 SR_50M.tif SR_50M-box.tif

And project into Mercator projection :

gdalwarp -overwrite -s_srs "+proj=longlat +datum=WGS84 +no_defs" -t_srs EPSG:3395 SR_50M-box.tif SR_50M-proj.tif

Finally I export to PNG :

gdal_translate -of PNG SR_50M-proj.tif SR_50M.png

Then I render everything using D3

   this.projection = d3.geoMercator()
        .translate([0, 0])
        .scale(1)


   this.rasterImage = new Image();
   this.rasterImage.src = raster;

   this.path = d3.geoPath(this.projection, this.ctx);
   this.bb = this.path.bounds(this.land);

   const s = 1 / Math.max((this.bb[1][0] - this.bb[0][0]) / this.getWidth(), (this.bb[1][1] - this.bb[0][1]) / this.getHeight());
    // transform
   const t = [(this.getWidth() - s * (this.bb[1][0] + this.bb[0][0])) / 2, (this.getHeight() - s * (this.bb[1][1] + this.bb[0][1])) / 2];

    // update projection
   this.projection
       .scale(s)
       .translate(t);

   this.raster_width = (this.bb[1][0] - this.bb[0][0]) * s;
   this.raster_height = (this.bb[1][1] - this.bb[0][1]) * s;

   this.rtranslate_x = (this.getWidth() - this.raster_width) / 2;
   this.rtranslate_y = (this.getHeight() - this.raster_height) / 2;

   this.ctx.beginPath();
   this.path(this.land);
   this.ctx.fill();

   this.ctx.save();
   this.ctx.globalAlpha = 0.8;
   this.ctx.translate(this.rtranslate_x, this.rtranslate_y);
   this.ctx.drawImage(this.rasterImage, 0, 0, this.raster_width, this.raster_height);
   this.ctx.restore();

However at the end, the geojson land and the raster are not aligned :

result

I tried to first project and then cut, or to use Pseudo-Mercator or Mercator when projection, but nothing works.. Anyone have an idea ?

Solution

  • A Mercator is usually clipped at roughly 85 degrees North/South (~85.05113 N/S) - as further than this you get a map that is taller than it is wide, and one that gets much much taller for every extra degree north/south included in the extent..

    D3 clips features using this limit:

    The spherical Mercator projection. Defines a default projection.clipExtent such that the world is projected to a square, clipped to approximately ±85° latitude.

    The northern bounds are fine, but the southern bounds of the geojson is -89.998926 degrees which you use to cut the image. But as D3 clips the geojson, your stretching the image by a different amount as compared with the geojson, hence the issue you see.

    The solution should be to clip the image to a bounds that is representative of the limits of what D3 will render for a Mercator (85.05113 degrees south) not the limits of the data itself.

    I haven't looked up how faithfully gdal implements EPSG:3395 as the definition provides for a projected bounds of 80 degrees south and 84 degrees north - though looking at the image, this doesn't appear to be an issue.

    You can also use the cleaner fitSize methods for D3 projections (d3v4+):

     projection.fitSize([width,height],geojsonObject)

    Which will set scale and translate for you given a provided width/height.