GeoTIFF Contours II Adapted from https://bl.ocks.org/mbostock/83c0be21dba7602ee14982b020b12f51 for https://github.com/d3/d3-geo-projection/issues/105 (See also https://beta.observablehq.com/@fil/netcdf)
draw = { const svg = DOM.svg(960, 500); var path = d3.geoPath(projection); d3.select(svg) .attr("stroke", "#000") .attr("stroke-width", 0.5) .attr("stroke-linejoin", "round") .selectAll("path") .data(data) .enter() .append("path") .attr("fill", function(d) { return color(d.value); }) .attr("d", path); return svg; }
data = { var values = rotate((await image.readRasters())[0]); color.domain(d3.extent(values)); var contours = d3.contours().size([n, m]); return contours(values).map(invert); // 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 > 0 ? { type: "Polygon", coordinates: p.coordinates, value: d.value } : { type: "Sphere", value: d.value }; }
load geotiff
tiff = GeoTIFF.fromUrl( "https://cors-anywhere.herokuapp.com/https://gist.githack.com/mbostock/83c0be21dba7602ee14982b020b12f51/raw/2300f3fa26bace2b1126e2ec096c1bc19ffc2024/sfctmp.tiff" )
image = tiff.getImage()
m = image.getHeight()
n = image.getWidth()
libraries
GeoTIFF = require("https://unpkg.com/geotiff@1.0.0-beta.6/dist/geotiff.bundle.min.js")
color = d3.scaleSequential(d3.interpolateMagma)
projection = d3.geoNaturalEarth().precision(0.1)
d3 = require("d3@5", "d3-geo-projection@2")