armollica/Makefile

## README.md

      
    Raw
  

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

  
## contours.json

      
Display the source blob

    
Display the rendered blob

    
    Raw
  

              contours.json
            
          
      Loading

      Sorry, something went wrong. Reload?
      Sorry, we cannot display this file.
      Sorry, this file is invalid so it cannot be displayed.
      
          Viewer requires iframe.
      
    
## index.html
<!DOCTYPE html>
<head>
<style>

/* See http://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

## thumbnail.png

      
    Raw
  

              thumbnail.png
	<!DOCTYPE html>
	<head>
	<style>

	/* See http://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>
	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