block by armollica 6559f04439cb8b0ccd30a69c3fc3976f

Maryland Contour Map

Full Screen

Recipe for making a contour map (see Makefile). Requires make, gdal and mapshaper.

index.html

<!DOCTYPE html>
<head>
<style>

/* See //soliton.vm.bytemark.co.uk/pub/cpt-city/wkp/template/tn/wiki-1.02.png.index.html */
.elevation-0 { fill: rgb(172,208,165); }
.elevation-1 { fill: rgb(148,191,139); }
.elevation-2 { fill: rgb(168,198,143); }
.elevation-3 { fill: rgb(189,204,150); }
.elevation-4 { fill: rgb(209,215,171); }
.elevation-5 { fill: rgb(225,228,181); }
.elevation-6 { fill: rgb(239,235,192); }
.elevation-7 { fill: rgb(232,225,182); }
.elevation-8 { fill: rgb(222,214,163); }
.elevation-9 { fill: rgb(211,202,157); }
.elevation-10 { fill: rgb(202,185,130); }
.elevation-11 { fill: rgb(195,167,107); }

</style>
</head>
<body>
<script src="https://d3js.org/d3.v4.min.js"></script>
<script src="https://d3js.org/topojson.v2.min.js"></script>
<script>

var width = 960,
		height = 500;

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

// EPSG:26985 (See https://github.com/veltman/d3-stateplane)
var projection = d3.geoConicConformal()
  .parallels([38 + 18 / 60, 39 + 27 / 60])
  .rotate([77, -37 - 40 / 60]);

var path = d3.geoPath()
	.projection(projection);

d3.json("contours.json", function(error, topology) {
	if (error) throw error;

	var contours = topojson.feature(topology, topology.objects.contours);

	projection.fitSize([width, height], contours);

	svg.selectAll("path").data(contours.features)
		.enter().append("path")
			.attr("class", function(d) { return "elevation-" + d.properties.DN; })
			.attr("d", path);
});

</script>
</body>
</html>

Makefile

GENERATED_FILES = \
	gz/statesp010g.shp_nt00938.tar.gz \
	shp/statesp010g.shp \
	shp/md-border.shp \
	zip/srtm_21_05.zip \
	tif/srtm_21_05.tif \
	tif/dem-reclassed.tif \
	tif/dem-sieved.tif \
	shp/contours-uncropped.shp \
	geojson/contours.json \
	contours.json

all: $(GENERATED_FILES)

# State borders
gz/statesp010g.shp_nt00938.tar.gz:
	mkdir -p $(dir $@)
	curl 'https://prd-tnm.s3.amazonaws.com/StagedProducts/Small-scale/data/Boundaries/statesp010g.shp_nt00938.tar.gz' -o $@.download
	mv $@.download $@

shp/statesp010g.shp: gz/statesp010g.shp_nt00938.tar.gz
	mkdir -p $(dir $@)
	tar -xzm -C $(dir $@) -f $<

# Maryland border
shp/md-border.shp: shp/statesp010g.shp
	mkdir -p $(dir $@)
	ogr2ogr -where 'STATE_ABBR = "MD"' $@ $<

# Raw DEM for big chunk of Mid-Atlantic, including Maryland
zip/srtm_21_05.zip:
	mkdir -p $(dir $@)
	curl 'http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/srtm_21_05.zip' -o $@.download
	mv $@.download $@

tif/srtm_21_05.tif: zip/srtm_21_05.zip
	mkdir -p $(dir $@)
	unzip -d $(dir $@) $<
	touch $@

# Bin pixels into contours
tif/dem-reclassed.tif: tif/srtm_21_05.tif
	mkdir -p $(dir $@)
	gdal_calc.py \
		-A tif/srtm_21_05.tif \
		--overwrite \
		--NoDataValue=-1 \
		--outfile=$@ \
		--calc "1  * logical_and(A >=   0, A <   50) + \
						2  * logical_and(A >=  50, A <  100) + \
						3  * logical_and(A >= 100, A <  200) + \
						4  * logical_and(A >= 200, A <  300) + \
						5  * logical_and(A >= 300, A <  400) + \
						6  * logical_and(A >= 400, A <  500) + \
						7  * logical_and(A >= 500, A <  600) + \
						8  * logical_and(A >= 600, A <  700) + \
						9  * logical_and(A >= 700, A <  800) + \
						10 * logical_and(A >= 800, A <  900) + \
						11 * logical_and(A >= 900, A < 1000)"

# Remove small raster polygons with sieve filter
tif/dem-sieved.tif: tif/dem-reclassed.tif
	mkdir -p $(dir $@)
	gdal_sieve.py -st 1000 "$<" $@

# Convert raster elevtion values into contour polygons (...takes forever)
shp/contours-uncropped.shp: tif/dem-sieved.tif
	mkdir -p $(dir $@)
	gdal_polygonize.py -f "ESRI Shapefile" "$<" $@

# Crop to Maryland border
# Simplify the geometry
# Fix "no data" issue (DN = DN === -1 ? 0 : DN)
# Dissolve each elevation group into single multi-polygon
geojson/contours.json: shp/contours-uncropped.shp shp/md-border.shp
	mkdir -p $(dir $@)
	mapshaper shp/contours-uncropped.shp \
		-clip shp/md-border.shp \
		-simplify 35% \
		-each 'DN = DN === -1 ? 0 : DN' \
		-dissolve DN \
		-o $@ format=geojson force

# Convert to TopoJSON
contours.json: geojson/contours.json
	mapshaper $< -o $@ format=topojson force

.PHONY: all