30:[{"id":"6954c8d84aa4a9b4de4dfcc4cbe86556571ecdcbbac8f587cc98207c9a613b03ad260859f796a61b2facf50b50038f113c8cd5f55c4fbe1e2c1bc785e8ac3805","href":"https://static.observableusercontent.com/files/6954c8d84aa4a9b4de4dfcc4cbe86556571ecdcbbac8f587cc98207c9a613b03ad260859f796a61b2facf50b50038f113c8cd5f55c4fbe1e2c1bc785e8ac3805","name":"sfctmp.tiff","type":"image/tiff","lastModified":1572392055885.025,"size":2083493}] 2e:{"title":"GeoTIFF contours II","cells":[{"id":0,"value":"

GeoTIFF contours II

D3 › Gallery

\n\n# GeoTIFF contours II\n\nThese contours are computed in equirectangular coordinates and then reprojected.","pinned":false,"hidden":false,"mode":"md"},{"id":31,"value":"chart = html``","pinned":true,"hidden":false,"mode":"ojs"},{"id":38,"value":"projection = d3.geoNaturalEarth1().precision(0.1)","pinned":true,"hidden":false,"mode":"ojs"},{"id":33,"value":"path = d3.geoPath(projection)","pinned":true,"hidden":false,"mode":"ojs"},{"id":29,"value":"contours = d3.contours().size([n, m])","pinned":true,"hidden":false,"mode":"ojs"},{"id":27,"value":"color = d3.scaleSequential(d3.extent(values), d3.interpolateMagma)","pinned":true,"hidden":false,"mode":"ojs"},{"id":16,"value":"image = tiff.getImage()","pinned":true,"hidden":false,"mode":"ojs"},{"id":15,"value":"values = rotate((await image.readRasters())[0])","pinned":true,"hidden":false,"mode":"ojs"},{"id":24,"value":"n = image.getWidth()","pinned":true,"hidden":false,"mode":"ojs"},{"id":22,"value":"m = image.getHeight()","pinned":true,"hidden":false,"mode":"ojs"},{"id":47,"value":"The *rotate* function rotates a GeoTIFF’s longitude from ${tex`[0, 360]`} to ${tex`[-180, +180]`}.","pinned":false,"hidden":false,"mode":"md"},{"id":42,"value":"function rotate(values) {\n var l = n >> 1;\n for (var j = 0, k = 0; j < m; ++j, k += n) {\n values.subarray(k, k + l).reverse();\n values.subarray(k + l, k + n).reverse();\n values.subarray(k, k + n).reverse();\n }\n return values;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":50,"value":"The *invert* function converts [*x*, *y*] in pixel coordinates to [*longitude*, *latitude*]. This assumes the source GeoTIFF is in equirectangular coordinates.\n\nInverting 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!","pinned":false,"hidden":false,"mode":"md"},{"id":45,"value":"function invert(d) {\n const shared = {};\n\n let p = {\n type: \"Polygon\",\n coordinates: d3.merge(d.coordinates.map(polygon => {\n return polygon.map(ring => {\n return ring.map(point => {\n return [point[0] / n * 360 - 180, 90 - point[1] / m * 180];\n }).reverse();\n });\n }))\n };\n\n // Record the y-intersections with the antimeridian.\n p.coordinates.forEach(ring => {\n ring.forEach(p => {\n if (p[0] === -180) shared[p[1]] |= 1;\n else if (p[0] === 180) shared[p[1]] |= 2;\n });\n });\n\n // Offset any unshared antimeridian points to prevent their stitching.\n p.coordinates.forEach(ring => {\n ring.forEach(p => {\n if ((p[0] === -180 || p[0] === 180) && shared[p[1]] !== 3) {\n p[0] = p[0] === -180 ? -179.9995 : 179.9995;\n }\n });\n });\n\n p = d3.geoStitch(p);\n\n // If the MultiPolygon is empty, treat it as the Sphere.\n return p.coordinates.length\n ? {type: \"Polygon\", coordinates: p.coordinates}\n : {type: \"Sphere\"};\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":13,"value":"tiff = FileAttachment(\"sfctmp.tiff\").arrayBuffer()\n .then(buffer => GeoTIFF.fromArrayBuffer(buffer))","pinned":true,"hidden":false,"mode":"ojs"},{"id":6,"value":"GeoTIFF = require(\"geotiff@2.0.7\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":5,"value":"d3 = require(\"d3@7\", \"d3-geo-projection@4\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":91,"value":"Or, using [Observable Plot](/plot/)’s concise API (see [complete example](https://observablehq.com/@observablehq/plot-geotiff-contours)):","pinned":false,"hidden":false,"mode":"md"},{"id":94,"value":"Plot.plot({\n projection: \"equal-earth\",\n color: { scheme: \"Magma\" },\n marks: [Plot.contour(values, {\n x: (_, i) => i % n / 2 - 180,\n y: (_, i) => 90 - Math.floor(i / n) / 2,\n fill: values,\n thresholds: 10,\n stroke: \"#000\",\n strokeWidth: 0.5,\n clip: \"sphere\"\n })]\n})","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W30","stdlib":"1","titling":"auto"}