11

Given a GIS raster with elevation data, How to design a topographic map in D3js ?

Is there any example of relief / topographic maps of cropped lands made using D3js ?


Not working: I explored the posibility of .tif > gdal_contour.py > .shp > topojson > d3js without success.

I use a makefile which contain all my commands. As my area of interest (France) is a crop of land areas, the gdal_contour.py approach generates broken isolines which does NOT create closed polygons. Also, the SVG end result fails. The only example of D3 topographic map I know is about Iceland, which, as an island, avoid this issue: cropping the country out of the world doesn't result in broken isolines.

enter image description here

nb: This project is part of the #Wikipedia #wikimaps project.

Hugolpz
  • 17,296
  • 26
  • 100
  • 187

2 Answers2

27

Topographic map now on D3js, with full makefile workflow ! See http://bl.ocks.org/hugolpz/6279966 (<= older code, compare to here on SO)

0. Requirements:

  • Geographic area: You may customize your geographic area of interest by editing one line within each of the 2 files : makefile#boxing and html#Geo-frame_borders with your own decimal coordinates fo W,N,E,S borders, something like:

    var WNES = { "target": "France", "W": -5.3, "N":51.6, "E": 10.2, "S": 41.0 };

  • Software: make, curl, unzip, gdal (include ogr, gdal_calc.py, gdal_polygonize.py), nodejs, topojson. Helpful: touch. The makefile then manage to download the sources, process them, and output a single topojson file that the D3js code provided can use.

1. Save into folder name: /topo_map/topo.mk

# topojsoning: 
final.json:  levels.json
    topojson --id-property none --simplify=0.5 -p name=elev -o final.json -- levels.json
    # simplification approach to explore further. Feedbacks welcome. 

# shp2jsoning:
levels.json: levels.shp
    ogr2ogr -f GeoJSON -where "elev < 10000" levels.json levels.shp

# merge
levels.shp: level0001.shp level0050.shp level0100.shp level0200.shp level0500.shp level1000.shp level2000.shp level3000.shp level4000.shp level5000.shp
    ogr2ogr levels.shp level0001.shp
    ogr2ogr -update -append levels.shp level0050.shp
    ogr2ogr -update -append levels.shp level0100.shp
    ogr2ogr -update -append levels.shp level0200.shp
    ogr2ogr -update -append levels.shp level0500.shp
    ogr2ogr -update -append levels.shp level1000.shp
    ogr2ogr -update -append levels.shp level2000.shp
    ogr2ogr -update -append levels.shp level3000.shp
    ogr2ogr -update -append levels.shp level4000.shp
    ogr2ogr -update -append levels.shp level5000.shp

# Polygonize slices:
level0001.shp: level0001.tif
    gdal_polygonize.py level0001.tif -f "ESRI Shapefile" level0001.shp level_0001 elev
level0050.shp: level0050.tif
    gdal_polygonize.py level0050.tif -f "ESRI Shapefile" level0050.shp level_0050 elev
level0100.shp: level0100.tif
    gdal_polygonize.py level0100.tif -f "ESRI Shapefile" level0100.shp level_0100 elev
level0200.shp: level0200.tif
    gdal_polygonize.py level0200.tif -f "ESRI Shapefile" level0200.shp level_0200 elev
level0500.shp: level0500.tif
    gdal_polygonize.py level0500.tif -f "ESRI Shapefile" level0500.shp level_0500 elev
level1000.shp: level1000.tif
    gdal_polygonize.py level1000.tif -f "ESRI Shapefile" level1000.shp level_1000 elev
level2000.shp: level2000.tif
    gdal_polygonize.py level2000.tif -f "ESRI Shapefile" level2000.shp level_2000 elev
level3000.shp: level3000.tif
    gdal_polygonize.py level3000.tif -f "ESRI Shapefile" level3000.shp level_3000 elev
level4000.shp: level4000.tif
    gdal_polygonize.py level4000.tif -f "ESRI Shapefile" level4000.shp level_4000 elev
level5000.shp: level5000.tif
    gdal_polygonize.py level5000.tif -f "ESRI Shapefile" level5000.shp level_5000 elev

# Raster slicing:
level0001.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0001.tif --calc="1*(A>0)"       --NoDataValue=0
level0050.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0050.tif --calc="50*(A>50)"      --NoDataValue=0
level0100.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0100.tif --calc="100*(A>100)"     --NoDataValue=0
level0200.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0200.tif --calc="200*(A>200)"     --NoDataValue=0
level0500.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0500.tif --calc="500*(A>500)"     --NoDataValue=0
level1000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level1000.tif --calc="1000*(A>1000)"     --NoDataValue=0
level2000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level2000.tif --calc="2000*(A>2000)"     --NoDataValue=0
level3000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level3000.tif --calc="3000*(A>3000)"     --NoDataValue=0
level4000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level4000.tif --calc="4000*(A>4000)"     --NoDataValue=0
level5000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level5000.tif --calc="5000*(A>5000)"     --NoDataValue=0

# boxing: 
crop.tif: ETOPO1_Ice_g_geotiff.tif
    gdal_translate -projwin -5.3 41.0 10.2 51.6 ETOPO1_Ice_g_geotiff.tif crop.tif
    # ulx uly lrx lry  // W S E N

# unzip:
ETOPO1_Ice_g_geotiff.tif: ETOPO1.zip
    unzip ETOPO1.zip
    touch ETOPO1_Ice_g_geotiff.tif

# download:
ETOPO1.zip:
    curl -o ETOPO1.zip 'http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/georeferenced_tiff/ETOPO1_Ice_g_geotiff.zip'

clean:
    rm `ls | grep -v 'zip' | grep -v 'Makefile'`
# Makefile v4b (@hugo_lz) 

2. Create data by running the makfile :

cd ./topo_map
make -f ./topo.mk

3. D3js & HTML code with auto-focus:

<!-- language: html -->
<style>
svg { border: 5px solid #333; background-color: #C6ECFF;}

/* TOPO */
path.Topo_1 { fill:#ACD0A5; stroke: #0978AB; stroke-width: 1px; }
path.Topo_50 {fill: #94BF8B; }
path.Topo_100 {fill: #BDCC96; }
path.Topo_200 {fill: #E1E4B5; }
path.Topo_500 {fill: #DED6A3; }
path.Topo_1000 {fill:#CAB982 ; }
path.Topo_2000 {fill: #AA8753; }
path.Topo_3000 {fill: #BAAE9A; }
path.Topo_4000 {fill: #E0DED8 ; }
path.Topo_5000 {fill: #FFFFFF ; }
.download { 
  background: #333; 
  color: #FFF; 
  font-weight: 900; 
  border: 2px solid #B10000; 
  padding: 4px; 
  margin:4px;
}
</style>
<body>
<script src="http://code.jquery.com/jquery-2.0.2.min.js"></script>
<script src="http://d3js.org/d3.v3.min.js"></script>
<script src="http://d3js.org/topojson.v1.min.js"></script>
<script>
// 1. -------------- SETTINGS ------------- //
// Geo-frame_borders in decimal ⁰: France
var WNES = { "W": -5.3, "N":51.6, "E": 10.2, "S": 41.0 };

// Geo values of interest :
var latCenter = (WNES.S + WNES.N)/2,
    lonCenter = (WNES.W + WNES.E)/2,
    geo_width = (WNES.E - WNES.W),
    geo_height= (WNES.N - WNES.S);
// HTML expected frame dimensions
var width  = 600,
    height = width * (geo_height / geo_width);

// Projection: projection, reset scale and translate
var projection = d3.geo.equirectangular()
      .scale(1)
      .translate([0, 0]);

// SVG injection:
var svg = d3.select("body").append("svg")
    .attr("width", width)
    .attr("height", height);

// Path
var path = d3.geo.path()
    .projection(projection)
    .pointRadius(4);

// Data (getJSON: TopoJSON)
d3.json("final.json", showData);

// 2. ---------- FUNCTION ------------- //
function showData(error, fra) {
    var Levels = topojson.feature(fra, fra.objects.levels);

// Focus area box compute for derive scale & translate.
// [​[left, bottom], [right, top]​] // E   W    N   S
var b = path.bounds(Levels),
    s = 1 / Math.max((b[1][0] - b[0][0]) / width, (b[1][1] - b[0][1]) / height),
    t = [(width - s * (b[1][0] + b[0][0])) / 2, (height - s * (b[1][1] + b[0][1])) / 2];

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

//Append Topo polygons
    svg.append("path")
        .datum(Levels)
        .attr("d", path)
    svg.selectAll(".levels")
        .data(topojson.feature(fra, fra.objects.levels).features)
      .enter().append("path")
        .attr("class", function(d) { return "Topo_" + d.properties.name; })
        .attr("data-elev", function(d) { return d.properties.name; })
        .attr("d", path)

}
</script>
<br />
<div>
    <a class="download ac-icon-download" href="javascript:javascript: (function () { var e = document.createElement('script'); if (window.location.protocol === 'https:') { e.setAttribute('src', 'https://raw.github.com/NYTimes/svg-crowbar/gh-pages/svg-crowbar.js'); } else { e.setAttribute('src', 'http://nytimes.github.com/svg-crowbar/svg-crowbar.js'); } e.setAttribute('class', 'svg-crowbar'); document.body.appendChild(e); })();"><!--⤋--><big>⇩</big> Download</a> -- Works on Chrome. Feedback me for others web browsers ?
</div>
<br />
</body>
</html>

4. Result will be precisely such: (applied to your area of interest)

enter image description here

If you publish map(s) online please share up the link :)

Note: encouraging +1 welcome.

Hugolpz
  • 17,296
  • 26
  • 100
  • 187
  • Requirements#Softwares: I installed these softwares a while ago, at different times. If someone get her/his hands on parts or all of the `sudo apt-get install` commands to run, thanks to share here or online :) – Hugolpz Aug 26 '13 at 06:35
  • I see the following error when I try the above Makefile: `Input file size is 21601, 10801 Computed -srcwin 10482 2940 930 -635 from projected window. Error: Computed -srcwin 10482 2940 930 -635 has negative width and/or height. *** [crop.tif] Error code 1` – Steve Wart Apr 18 '14 at 18:23
  • The task `crop.tif:` is failling. Seems your width and/or height is negative. Also, check your `-projwin` values and their places within the line `gdal_translate -projwin -5.3 41.0 10.2 51.6 ETOPO1_Ice_g_geotiff.tif crop.tif`. The order is West, North, East, South borders, but may need a change if you map the pacific and cross the [180th meridian](https://en.wikipedia.org/wiki/180th_meridian). – Hugolpz Apr 19 '14 at 11:12
  • Thanks. I am using the Makefile exactly as it appears above. I'll fiddle with the coordinates and see if I can figure out what's wrong – Steve Wart Apr 19 '14 at 15:03
  • @SteveWart : did you got something nice? – Hugolpz Feb 01 '15 at 18:04
  • Here is a link http://tardis.mandeluna.com/topomap/ I seem to recall my problem was related to the WNES parameters. Eventually I used the following { "W": 9.568, "N":41.99, "E": 19.73, "S": 36.34 }; Thanks for your help. I've uploaded the final files if anyone is interested in the approach, but not any of the rather large data files. My Makefile is based on the one you posted http://tardis.mandeluna.com/topomap/Makefile - I last used this code in May 2014 and it probably has errors. Use at your own risk. – Steve Wart Feb 03 '15 at 04:59
  • Cool. The edge are crispy, cool get better if using higher quality topographic data (SRTM). – Hugolpz Feb 03 '15 at 09:22
  • Awesome code; worked great for me. However, for small areas (map of a city or region), the NGDC/NOAA TIFF which you use in the above results in very pixelated features. Do you have any experience or suggestions in adapting this to work with other elevation data sources (SRTM/DEM)? I'm working on maps of San Francisco and Manhattan :) – emunsing Aug 04 '15 at 06:31
  • 90m resolution http://dwtkns.com/srtm/ , if you want general shapes for a 20*20km area. For smaller (manathan island, cities), you will need 31m or lower. The USA has such public data, it should take you few mins to find something. – Hugolpz Aug 04 '15 at 08:58
  • Will this method work with a global map or are there limitations on the size of the area? – AnnanFay Sep 12 '18 at 20:39
  • 1
    @Annan: Yes, should work. 1) You could remove the task #boxing since you keep it all, and edit the rest of the file accordingly. If so, you will process a HUGE TIFF FILE. So I would also suggest you to do a simple resize over ETOPO1_Ice_g_geotiff.tif so it gets about 2~4000px width. Then you can play. – Hugolpz Sep 17 '18 at 09:34
1

If anyone is looking for an update, here's the build code I got running as of today. Required me to manually download the .zip file and move it into the topo_map directory, and then a few changes (noted in bold):

# topojsoning (USE GEO2TOPO not TOPOJSON): 
final.json: levels.json
    geo2topo --id-property none --simplify=0.5 -p name=elev -o final.json -- levels.json
    # simplification approach to explore further. Feedbacks welcome. 

# shp2jsoning:
levels.json: levels.shp
    ogr2ogr -f GeoJSON -where "elev < 10000" levels.json levels.shp

# merge
levels.shp: level0001.shp level0050.shp level0100.shp level0200.shp level0500.shp level1000.shp level2000.shp level3000.shp level4000.shp level5000.shp
    ogr2ogr levels.shp level0001.shp
    ogr2ogr -update -append levels.shp level0050.shp
    ogr2ogr -update -append levels.shp level0100.shp
    ogr2ogr -update -append levels.shp level0200.shp
    ogr2ogr -update -append levels.shp level0500.shp
    ogr2ogr -update -append levels.shp level1000.shp
    ogr2ogr -update -append levels.shp level2000.shp
    ogr2ogr -update -append levels.shp level3000.shp
    ogr2ogr -update -append levels.shp level4000.shp
    ogr2ogr -update -append levels.shp level5000.shp

# Polygonize slices:
level0001.shp: level0001.tif
    gdal_polygonize.py level0001.tif -f "ESRI Shapefile" level0001.shp level_0001 elev
level0050.shp: level0050.tif
    gdal_polygonize.py level0050.tif -f "ESRI Shapefile" level0050.shp level_0050 elev
level0100.shp: level0100.tif
    gdal_polygonize.py level0100.tif -f "ESRI Shapefile" level0100.shp level_0100 elev
level0200.shp: level0200.tif
    gdal_polygonize.py level0200.tif -f "ESRI Shapefile" level0200.shp level_0200 elev
level0500.shp: level0500.tif
    gdal_polygonize.py level0500.tif -f "ESRI Shapefile" level0500.shp level_0500 elev
level1000.shp: level1000.tif
    gdal_polygonize.py level1000.tif -f "ESRI Shapefile" level1000.shp level_1000 elev
level2000.shp: level2000.tif
    gdal_polygonize.py level2000.tif -f "ESRI Shapefile" level2000.shp level_2000 elev
level3000.shp: level3000.tif
    gdal_polygonize.py level3000.tif -f "ESRI Shapefile" level3000.shp level_3000 elev
level4000.shp: level4000.tif
    gdal_polygonize.py level4000.tif -f "ESRI Shapefile" level4000.shp level_4000 elev
level5000.shp: level5000.tif
    gdal_polygonize.py level5000.tif -f "ESRI Shapefile" level5000.shp level_5000 elev

# Raster slicing:
level0001.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0001.tif --calc="1*(A>0)"       --NoDataValue=0
level0050.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0050.tif --calc="50*(A>50)"      --NoDataValue=0
level0100.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0100.tif --calc="100*(A>100)"     --NoDataValue=0
level0200.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0200.tif --calc="200*(A>200)"     --NoDataValue=0
level0500.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level0500.tif --calc="500*(A>500)"     --NoDataValue=0
level1000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level1000.tif --calc="1000*(A>1000)"     --NoDataValue=0
level2000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level2000.tif --calc="2000*(A>2000)"     --NoDataValue=0
level3000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level3000.tif --calc="3000*(A>3000)"     --NoDataValue=0
level4000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level4000.tif --calc="4000*(A>4000)"     --NoDataValue=0
level5000.tif: crop.tif
    gdal_calc.py -A crop.tif --outfile=level5000.tif --calc="5000*(A>5000)"     --NoDataValue=0

# boxing: 
crop.tif: ETOPO1_Ice_g_geotiff.tif
    gdal_translate -projwin -84.9 47.0 -69.9 33.7 ETOPO1_Ice_g_geotiff.tif crop.tif
    # ulx uly lrx lry  // W N E S <- Coordinate order
# unzip:
ETOPO1_Ice_g_geotiff.tif: ETOPO1.zip
    unzip ETOPO1.zip
    touch ETOPO1_Ice_g_geotiff.tif

# download:
#ETOPO1.zip:
 #   curl -o ETOPO1.zip 'http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/georeferenced_tiff/ETOPO1_Ice_g_geotiff.zip'

clean:
    rm `ls | grep -v 'zip' | grep -v 'Makefile'`
# Makefile v4b (@Lopez_lz) 
  • Thanks Nick. The differences are visible there https://www.diffchecker.com/lFTXbr1A . Could you clarify why are these edits in your answer ? As for the **CLI command** ```topojson``` => ```geo2json```, was there a CLI change ? For the **download url**, ideally, whenever the url provided fails you just have to find and update the url within the makefile. For the ```WSEN``` => ```WNES```, the [gdal_translate manual](http://gdal.org/gdal_translate.html) seems to confirm my initial configuration. Yours may works as well since gdal is smart :) – Hugolpz Dec 21 '16 at 11:09
  • Yeah I'll edit my response to reflect the changes. There are two main changes are topojson => geo2json and the coordinate order for cropping. – Nick Malawskey Dec 23 '16 at 14:54