Create a gist now

Instantly share code, notes, and snippets.

@mjhoy /Makefile
Last active Feb 16, 2018

Embed
What would you like to do?
Costa Rica shaded relief

Shaded Relief

This is an example of how to create a shaded relief raster with a vector map overlay (using SVG and d3.js).

Step 1 was to create the raster. I used tiled GeoTiffs from the SRTM project, downloading four tiles that completed a map of Costa Rica. To combine the tiff files into a single raster with the correct projection and dimensions, I used gdalwarp:

gdalwarp \
  -r lanczos \
  -te -250000 -156250 250000 156250 \
  -t_srs "+proj=aea +lat_1=8 +lat_2=11.5 +lat_0=9.7 +lon_0=-84.2 +x_0=0 +y_0=0" \
  -ts 960 0 \
  srtm_19_10.tif \
  srtm_20_10.tif \
  srtm_19_11.tif \
  srtm_20_11.tif \
  relief.tiff

The t_srs option sets an albers equal area projection that will center on Costa Rica. The te option defines the extent of the map, using SRS coordinates. I don't fully understand how this works and used some trial and error. Note that the x/y has a ratio of 1.6, the same as the intended output resolution (960x600).

Note that the projection here mirrors the projection set in index.html.

Step 2 is to create, from this GeoTiff file, two images: one, grayscale, that represents "shade" — given a certain direction of sunlight, it simulates the effect of light on the relief map:

gdaldem \
  hillshade \
  relief.tiff \
  hill-relief-shaded.tiff \
  -z 4 -az 20

The second image is a "color relief" that maps certain colors to certain elevations. The color_relief.txt file provides this information in the format: elevation r g b.

gdaldem \
  color-relief \
  relief.tiff \
  color_relief.txt \
  hill-relief-c.tiff \

These files are combined using the program hsv_merge.py:

hsv_merge.py \
  hill-relief-c.tiff \
  hill-relief-shaded.tiff \
  hill-relief-merged.tiff
65535 255 255 255
5800 254 254 254
3000 121 117 10
1500 151 106 47
800 127 166 122
500 213 213 149
1 201 213 166
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
<!DOCTYPE html>
<meta charset='utf-8'>
<style>
/* Style here... */
path#CRI {
fill: none;
stroke: #000;
}
image.bg {
opacity: 0.2;
}
</style>
<div></div>
<script src="http://d3js.org/d3.v3.min.js"></script>
<script src="http://d3js.org/topojson.v0.min.js"></script>
<script>
var width = 960,
height = 600;
var projection = d3.geo.albers()
.center([0, 9.7])
.rotate([84.2,0])
.parallels([8, 11.5])
//.scale(10200)
.scale(12240)
.translate([width / 2, height / 2]);
var path = d3.geo.path()
.projection(projection);
var svg = d3.select("body").append("svg")
.attr("class", "map")
.attr("width", width)
.attr("height", height);
d3.json("costarica_min_topo.json", function(error, data) {
var costarica = topojson.object(data, data.objects.costarica);
svg.append("image")
//.attr("clip-path", "url(#clip)")
.attr("xlink:href", "hill-relief.jpg")
.attr("width", width)
.attr("height", height)
.attr("class", "bg");
svg.selectAll(".cr-subunit")
.data(costarica.geometries)
.enter().append("path")
.attr("id", function(d) { return "CRI"; })
.attr("d", path);
svg.append("clipPath")
.attr("id", "clip")
.append("use")
.attr("xlink:href", "#CRI");
svg.append("image")
.attr("clip-path", "url(#clip)")
.attr("xlink:href", "hill-relief.jpg")
.attr("width", width)
.attr("height", height)
.attr("class", "fg");
});
// Script here...
</script>
all: hill-relief.jpg costarica_min_topo.json
# -------------
# Relief raster
# -------------
#
# Notice the `zip` file requirements here have no download.
# You will need to search for them online. They are from the
# SRTM project: http://www2.jpl.nasa.gov/srtm/
# (which appears to have multiple versions of files).
srtm_20_11.tif: srtm_20_11.zip
unzip srtm_20_11.zip
rm readme.txt
touch srtm_20_11.tif
srtm_20_10.tif: srtm_20_10.zip
unzip srtm_20_10.zip
rm readme.txt
touch srtm_20_10.tif
srtm_19_11.tif: srtm_19_11.zip
unzip srtm_19_11.zip
rm readme.txt
touch srtm_19_11.tif
srtm_19_10.tif: srtm_19_10.zip
unzip srtm_19_10.zip
rm readme.txt
touch srtm_19_10.tif
# Generate, from the four tiles, a relief GeoTiff file with the correct
# projection. Right now I'm finding the correct `-te` values by trial and error.
# -te -300000 -250000 300000 250000
relief.tiff: srtm_19_10.tif srtm_20_10.tif srtm_20_11.tif srtm_19_11.tif
gdalwarp \
-r lanczos \
-te -250000 -156250 250000 156250 \
-t_srs "+proj=aea +lat_1=8 +lat_2=11.5 +lat_0=9.7 +lon_0=-84.2 +x_0=0 +y_0=0" \
-ts 960 0 \
srtm_19_10.tif \
srtm_20_10.tif \
srtm_19_11.tif \
srtm_20_11.tif \
relief.tiff
# Shaded relief.
hill-relief-shaded.tiff: relief.tiff
gdaldem \
hillshade \
relief.tiff \
hill-relief-shaded.tiff \
-z 4 -az 20
# Colorized elevation.
hill-relief-c.tiff: relief.tiff
gdaldem \
color-relief \
relief.tiff \
color_relief.txt \
hill-relief-c.tiff \
# Merge the colorized and shaded rasters. Requires `hsv_merge.py`:
# http://fwarmerdam.blogspot.com/2010/01/hsvmergepy.html
hill-relief-merged.tiff: hill-relief-c.tiff hill-relief-shaded.tiff
~/bin/hsv_merge.py \
hill-relief-c.tiff \
hill-relief-shaded.tiff \
hill-relief-merged.tiff
hill-relief.jpg: hill-relief-merged.tiff
convert hill-relief-merged.tiff hill-relief.jpg
# --------------------
# Costa Rica Geography
# --------------------
# Data from http://gadm.org/
CRI_adm.zip:
curl -o CRI_adm.zip http://gadm.org/data/shp/CRI_adm.zip
CRI_adm0.shp: CRI_adm.zip
unzip CRI_adm.zip
touch CRI_adm0.shp
costarica.json: CRI_adm0.shp
ogr2ogr \
-f GeoJSON \
costarica.json \
CRI_adm0.shp
# Require topojson: https://github.com/mbostock/topojson
# (this minifies/simplifies the data)
costarica_min_topo.json: costarica.json
topojson \
-p name=NAME \
-p name \
-q 1e4 \
-o costarica_min_topo.json \
costarica.json
.PHONY: clean
# Clean does not remove downloaded files.
clean:
rm -f hill-relief*
rm -f relief*
rm -f srtm_19_10.hdr
rm -f srtm_19_10.t*
rm -f srtm_19_11.hdr
rm -f srtm_19_11.t*
rm -f srtm_20_10.hdr
rm -f srtm_20_10.t*
rm -f srtm_20_11.hdr
rm -f srtm_20_11.t*
rm -f CRI_adm0*
rm -f CRI_adm1*
rm -f CRI_adm2*
rm -f CRI_readme.txt
rm -f costarica.json
rm -f costarica_min_topo.json
@mjhoy

This comment has been minimized.

Show comment
Hide comment
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment