armollica · September 9, 2021 13:25
diff --git a/.block b/.block
 height: 1000
diff --git a/README.md b/README.md
diff --git a/index.html b/index.html
 <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>
diff --git a/make-borders.sh b/make-borders.sh
 #!/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
diff --git a/make-shaded-relief.sh b/make-shaded-relief.sh
 #!/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
diff --git a/shaded-relief.png b/shaded-relief.png
diff --git a/thumbnail.png b/thumbnail.png
diff --git a/wi.json b/wi.json
	<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