Skip to content

Instantly share code, notes, and snippets.

@armollica
Last active September 9, 2021 13:25
Show Gist options
  • Save armollica/47357a627782fe4d6ca447d9bfdabf1d to your computer and use it in GitHub Desktop.
Save armollica/47357a627782fe4d6ca447d9bfdabf1d to your computer and use it in GitHub Desktop.
Shaded Relief Map
height: 1000

Shaded relief map of Wisconsin.

The recipe for creating the shaded relief raster image came largely from Thomas Thoren's great blog post Making a New York Times Map. Derek Watkin's GDAL cheat sheet and SRTM Tile Grabber were also very useful.

This recipe is in the two bash scripts in this gist repo. First, make-borders.sh creates the Wisconsin state and county borders vector data. Second, make-shaded-relief.sh creates the shaded relief raster data from NASA's STRM data. The scripts need to be run in that order. The borders data is used to clip the raster data.

<html>
<head>
<style>
label {
position: absolute;
top: 15px;
right: 15px;
font-family: sans-serif;
font-size: 12px;
}
.state-interior {
fill: #fff;
}
.state-border {
fill: none;
stroke: lightslategrey;
}
.raster {
pointer-events: none;
fill: none;
}
.charcoalified .state-interior {
fill: #000;
}
.charcoalified .state-border {
stroke: none;
}
</style>
</head>
<body>
<label>
<input type="checkbox">
Charcoalify
</label>
<script src="https://d3js.org/d3.v3.min.js" charset="utf-8"></script>
<script src="//d3js.org/topojson.v1.min.js"></script>
<script>
var width = 960,
height = 1000;
// Path generator that works with preprojected geo data
// (see http://bl.ocks.org/mbostock/5663666)
var path = d3.geo.path()
.projection(matrix(1, 0, 0, -1, 0, height))
var svg = d3.select("body").append("svg")
.attr("width", width)
.attr("height", height);
d3.json("wi.json", function(error, wi) {
if (error) throw error;
var counties = topojson.feature(wi, wi.objects.counties),
state = topojson.merge(wi, wi.objects.counties.geometries);
// Fit the geography to the canvas (i.e., scale and translate)
var b = path.bounds(counties),
s = .975 / 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];
path.projection(matrix(s, 0, 0, -s, t[0], t[1]));
// Draw state interior
svg.append("path").datum(state)
.attr("class", "state-interior")
.attr("d", path);
// Draw shaded relief from raster data
var raster_width = (b[1][0] - b[0][0]) * s;
var raster_height = (b[1][1] - b[0][1]) * s;
var rtranslate_x = (width - raster_width) / 2;
var rtranslate_y = (height - raster_height) / 2;
svg.append("image")
.attr("xlink:href", "shaded-relief.png")
.attr("class", "raster")
.attr("width", raster_width)
.attr("height", raster_height)
.attr("transform",
"translate(" + rtranslate_x + ", " + rtranslate_y + ")");
// Draw state border from vector data
// (drawn after raster so it overlays)
svg.append("path").datum(state)
.attr("class", "state-border")
.attr("d", path);
// TODO: The raster data and vector data don't match up perfectly here
// This may be due to the fact that this state border is from merged
// county data and the raster image was clipped to a state border from
// a different vector file. May also be due to the projection issues.
});
// Charcoalify checkbox
var charcoalified = false;
d3.select("input")
.on("change", function() {
svg.classed("charcoalified", charcoalified = !charcoalified);
});
function matrix(a, b, c, d, tx, ty) {
return d3.geo.transform({
point: function(x, y) { this.stream.point(a * x + b * y + tx, c * x + d * y + ty); }
});
}
</script>
</body>
</html>
#!/usr/bin/env bash
#_______________________________________________________________________________
# Download and unzip shapefiles of US states and counties from Census
mkdir -p borders/zips
curl \
-o borders/zips/states.zip \
"http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_5m.zip"
curl \
-o borders/zips/counties.zip \
"ftp://ftp2.census.gov/geo/tiger/TIGER2015/COUNTY/tl_2015_us_county.zip"
unzip borders/zips/states.zip -d borders/us-states
unzip borders/zips/counties.zip -d borders/us-counties
#_______________________________________________________________________________
# Extract Wisconsin state and counties
# Filter states to Wisconsin
mkdir borders/wi-state
ogr2ogr \
-where "STATEFP = '55'" \
borders/wi-state/wi-state.shp \
borders/us-states/cb_2015_us_state_5m.shp
# Clip to WI state boundaries, filter to WI counties
mkdir borders/wi-counties
ogr2ogr \
-clipsrc borders/wi-state/wi-state.shp \
-where "STATEFP = '55'" \
borders/wi-counties/wi-counties.shp \
borders/us-counties/tl_2015_us_county.shp
#_______________________________________________________________________________
# Make Wisconsin counties ready for web display (counties.json)
# Convert to GeoJSON, reproject to latest Wisconsin Transverse Mercator and
# simplify the geometry a bit to shrink the filesize
ogr2ogr \
-f "GeoJSON" \
-t_srs "EPSG:3701" \
borders/counties.json \
borders/wi-counties/wi-counties.shp
# Convert to Topojson
topojson \
-o wi.json \
-- borders/counties.json
#!/usr/bin/env bash
#_______________________________________________________________________________
# Download elevation data
# URLs for the zip files containing the DEM GeoTiffs. Data are from NASA's
# Shuttle Radar Typography Mission (SRTM). Derek Watkin's SRTM Tile Grabber is a
# good resource for finding which files you need
urls="http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_20_04.zip"
urls="${urls} http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_20_03.zip"
urls="${urls} http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_19_04.zip"
urls="${urls} http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_19_03.zip"
urls="${urls} http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_18_03.zip"
urls="${urls} http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_18_04.zip"
mkdir -p dem/zips
cd dem/zips && { curl --remote-name-all $urls ; cd -; }
#_______________________________________________________________________________
# Unzip all GeoTiffs
unzip 'dem/zips/*.zip' \
-x readme.txt \
-d dem/unmerged
#_______________________________________________________________________________
# Merge data into a single GeoTiff
gdal_merge.py \
-o dem/merged.tif \
-init "255" \
dem/unmerged/*.tif
#_______________________________________________________________________________
# Reproject to latest Wisconsin Transverse Mercator
# (see http://spatialreference.org/ref/epsg/3701/)
gdalwarp \
-t_srs "EPSG:3701" \
-r bilinear \
dem/merged.tif \
dem/reprojected.tif
#_______________________________________________________________________________
# Create shaded relief
# Create raw hillshade greyscale
gdaldem hillshade \
dem/reprojected.tif \
dem/hillshade-raw.tif \
-s 2 \
-z 5 \
-alt 45 \
-compute_edges
# Crop to state border
gdalwarp \
-cutline borders/wi-state/wi-state.shp \
-crop_to_cutline \
-dstalpha \
dem/hillshade-raw.tif \
dem/hillshade.tif
## Fine tune hillshade
# Flat areas should be white (greyscale value = 255) and
# transparent (opacity value = 0). Need to find the cutoff
# value for flat land. In this case that value is 181. To
# find cutoff value use this command:
# gdallocationinfo -xml -wgs84 input.tif <lon> <lat>
# where the lon/lat point to a flat area. If
# the raster data was clipped then the area outside the clip
# is "flat" and is a good place to test.
# Create greyscale band
gdal_calc.py \
-A dem/hillshade.tif \
--outfile=dem/hillshade-greyscale.tif \
--calc="255*(A>180) + A*(A<=180)"
# Create opacity band
gdal_calc.py \
-A dem/hillshade.tif \
--outfile=dem/hillshade-opacity.tif \
--calc="0*(A>180) + (256-A)*(A<=180)"
# Combine into VRT as two separate bands
gdalbuildvrt \
-separate dem/hillshade.vrt \
dem/hillshade-greyscale.tif \
dem/hillshade-opacity.tif
# Convert this VRT back to TIF combining the greyscale
# and opacity bands
gdal_translate \
-co COMPRESS=LZW \
-co ALPHA=YES \
dem/hillshade.vrt \
dem/shaded-relief.tif
## Convert TIF to PNG
convert \
-resize x1200 \
dem/shaded-relief.tif \
shaded-relief.png
View raw

(Sorry about that, but we can’t show files that are this big right now.)

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