Skip to content

Instantly share code, notes, and snippets.

@mbostock

mbostock/.block

Last active November 19, 2023 22:53
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save mbostock/83c0be21dba7602ee14982b020b12f51 to your computer and use it in GitHub Desktop.
Save mbostock/83c0be21dba7602ee14982b020b12f51 to your computer and use it in GitHub Desktop.
Download ZIP
GeoTIFF Contours II
Raw
.block
license: gpl-3.0
border: no
redirect: https://observablehq.com/@d3/geotiff-contours-ii
Raw
README.md
Raw
index.html
<!DOCTYPE html>
<svg width="960" height="500"></svg>
<script src="https://unpkg.com/d3-array@1"></script>
<script src="https://unpkg.com/d3-contour@1"></script>
<script src="https://unpkg.com/d3-collection@1"></script>
<script src="https://unpkg.com/d3-color@1"></script>
<script src="https://unpkg.com/d3-dispatch@1"></script>
<script src="https://unpkg.com/d3-geo@1"></script>
<script src="https://unpkg.com/d3-geo-projection@2"></script>
<script src="https://unpkg.com/d3-interpolate@1"></script>
<script src="https://unpkg.com/d3-request@1"></script>
<script src="https://unpkg.com/d3-selection@1"></script>
<script src="https://unpkg.com/d3-scale@1"></script>
<script src="https://unpkg.com/geotiff@0.4/dist/geotiff.browserify.min.js"></script>
<script>
d3.request("sfctmp.tiff").responseType("arraybuffer").get(function(error, request) {
if (error) throw error;
var tiff = GeoTIFF.parse(request.response),
image = tiff.getImage(),
m = image.getHeight(),
n = image.getWidth(),
values = rotate(image.readRasters()[0]);
var color = d3.scaleSequential(d3.interpolateMagma)
.domain(d3.extent(values));
var projection = d3.geoNaturalEarth()
.precision(0.1);
var path = d3.geoPath(projection);
var contours = d3.contours()
.size([n, m]);
d3.select("svg")
.attr("stroke", "#000")
.attr("stroke-width", 0.5)
.attr("stroke-linejoin", "round")
.selectAll("path")
.data(contours(values).map(invert))
.enter().append("path")
.attr("fill", function(d) { return color(d.value); })
.attr("d", path);
// Rotate a GeoTIFF’s longitude from [0, 360] to [-180, +180].
function rotate(values) {
var l = n >> 1;
for (var j = 0, k = 0; j < m; ++j, k += n) {
values.subarray(k, k + l).reverse();
values.subarray(k + l, k + n).reverse();
values.subarray(k, k + n).reverse();
}
return values;
}
// Invert the pixel coordinates to [longitude, latitude]. This assumes the
// source GeoTIFF is in equirectangular coordinates.
//
// Note that inverting the projection breaks the polygon ring associations:
// holes are no longer inside their exterior rings. Fortunately, since the
// winding order of the rings is consistent and we’re now in spherical
// coordinates, we can just merge everything into a single polygon!
function invert(d) {
var shared = {};
var p = {
type: "Polygon",
coordinates: d3.merge(d.coordinates.map(function(polygon) {
return polygon.map(function(ring) {
return ring.map(function(point) {
return [point[0] / n * 360 - 180, 90 - point[1] / m * 180];
}).reverse();
});
}))
};
// Record the y-intersections with the antimeridian.
p.coordinates.forEach(function(ring) {
ring.forEach(function(p) {
if (p[0] === -180) shared[p[1]] |= 1;
else if (p[0] === 180) shared[p[1]] |= 2;
});
});
// Offset any unshared antimeridian points to prevent their stitching.
p.coordinates.forEach(function(ring) {
ring.forEach(function(p) {
if ((p[0] === -180 || p[0] === 180) && shared[p[1]] !== 3) {
p[0] = p[0] === -180 ? -179.9995 : 179.9995;
}
});
});
p = d3.geoStitch(p);
// If the MultiPolygon is empty, treat it as the Sphere.
return p.coordinates.length
? {type: "Polygon", coordinates: p.coordinates, value: d.value}
: {type: "Sphere", value: d.value};
}
});
</script>
Raw
sfctmp.tiff
This file has been truncated, but you can view the full file.
View raw
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment