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
Owner

mjhoy commented Apr 3, 2013

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