31:T510,{ if(map == "3DEP slope pct png"){ return d3.image(wmsSlopeUrlPng); } if(map == "3DEP slope pct geotiff"){ return plottyRenderer(slopeValues, rangeMin, rangeMax, colorramp); } if(map == "3DEP slope deg geotiff"){ return plottyRenderer(slopeDegValues, rangeMin, rangeMax, colorramp); } if(map == "3DEP aspect deg geotiff"){ return plottyRenderer(aspectDegValues, rangeMin, rangeMax, colorramp); } if(map == "3DEP aspect-slope png"){ return d3.image(wmsAspectUrlPng); } if(map == "3DEP aspect-slope geotiff"){ return plottyRenderer(aspectValues, rangeMin, rangeMax, colorramp); } if(map == "NLCD land cover png"){ return d3.image(wmsLandCoverUrlPng); } if(map == "NLCD land cover geotiff"){ return plottyRenderer(landCoverValues, rangeMin, rangeMax, colorramp); } if(map == "3DEP elevation-hillshade png"){ return d3.image(wmsElevationHillshadeUrlPng); } if(map == "3DEP elevation-hillshade geotiff"){ return plottyRenderer(elevHillshadeValues, rangeMin, rangeMax, colorramp); } if(map == "3DEP elevation geotiff"){ return plottyRenderer(elevationValues, rangeMin, rangeMax, colorramp); } if(map == "Mapzen elevation geotiff"){ return plottyRenderer(mapzenValues, rangeMin, rangeMax, colorramp);; } }32:T61f,viewof tileCoords = { const coords = { x: 28, y: 52, z: 7 }; const pan = html``; pan.up.onclick = () => move(0, 0, -1); pan.left.onclick = () => move(0, -1, 0); pan.right.onclick = () => move(0, 1, 0); pan.down.onclick = () => move(0, 0, 1); pan.zoomIn.onclick = () => move(1, 0, 0); pan.zoomOut.onclick = () => move(-1, 0, 0); function move(dz, dx, dy) { coords.x += dx; coords.y += dy; if (dz > 0 && coords.z < 14) { coords.z += dz; coords.x *= 2; coords.y *= 2; } else if (dz < 0 && coords.z > 4) { coords.z += dz; coords.x = Math.floor(coords.x / 2); coords.y = Math.floor(coords.y / 2); } pan.dispatchEvent(new CustomEvent("input")); } pan.value = coords; return pan; }33:T73b,html``34:[] 2f:{"title":"USGS 3DEP maps with adjustable colors","cells":[{"id":0,"value":"md`# USGS 3DEP maps with adjustable colors`","pinned":false,"hidden":false,"mode":"ojs"},{"id":927,"value":"md`Toggle through these high-res elevation-derived map tiles from the [USGS 3DEP program](https://www.usgs.gov/core-science-systems/ngp/3dep). Compare the png versions against the color-adjustable geotiff versions.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":277,"value":"viewof map = select({title:\"Select Map:\",\n options:[\"3DEP slope pct geotiff\", \"3DEP slope deg geotiff\", \"3DEP aspect deg geotiff\", \"3DEP aspect-slope geotiff\", \"3DEP elevation-hillshade geotiff\", \"3DEP elevation geotiff\", \"Mapzen elevation geotiff\",\"3DEP slope pct png\", \"3DEP aspect-slope png\", \"3DEP elevation-hillshade png\", \"NLCD land cover png\", \"NLCD land cover geotiff\"]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":641,"value":"viewof colorramp = select({title: \"Color Ramp (for geotiffs):\",\n options: [\"viridis\",\t\"inferno\",\t\"turbo\", \"rainbow\", \"jet\",\t\"hsv\", \"hot\",\t\"cool\",\t\"spring\", \"summer\",\t\"autumn\",\t\"winter\", \"bone\", \"copper\", \"greys\", \"yignbu\",\t\"greens\",\t\"yiorrd\", \"bluered\",\t\"rdbu\",\t\"picnic\", \"portland\",\t\"blackbody\",\t\"earth\", \"electric\",\t\"magma\",\t\"plasma\"]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":879,"value":"viewof rangeMin = number({value: range[0], title: \"Min color range (for geotiffs)\", submit: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":887,"value":"viewof rangeMax = number({value: range[1], title: \"Max color range (for geotiffs)\", submit: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":284,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":426,"value":"\"Z:\"+tileCoords.z+\", X:\"+tileCoords.x+\", Y:\"+tileCoords.y","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":107,"value":"md`These cool CSS navigation buttons are borrowed from [Jeshurun's notebook](https://observablehq.com/@jjhembd/vector-tile-rendering) and the handy selectors are from [Jeremy Ashkenas](https://observablehq.com/@jashkenas/inputs)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":91,"value":"md` [3DEP LIDAR](https://www.usgs.gov/core-science-systems/ngp/3dep) elevation data can be accessed from the [National Map](https://www.usgs.gov/core-science-systems/national-geospatial-program/national-map) using a [WMS service](https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetCapabilities&service=WMS) which returns a map within a bounding box. To use this service with other tiled maps, we need to be able to request a tile.

This notebook converts a tile z/x/y to a bounding box and requests data from the wms service. Geotiffs are read using [geotiffjs](https://github.com/geotiffjs/geotiff.js/) and then rendered using [plotty](http://santilland.github.io/plotty/index.html). The nice thing about using geotiffs is that we can adjust the colorbar and range.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":94,"value":"md`### Step1:\nConvert tile z/x/y to a bounding box. Equations are from the [OpenStreetMap Wiki](https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_bounding_box)` ","pinned":false,"hidden":false,"mode":"ojs"},{"id":14,"value":"bb=tiletobbox(tileCoords.z,tileCoords.x,tileCoords.y)","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"function tiletobbox(z,x,y) {\n \n var maxy = latToY(tile2lat(y, z));\n var miny = latToY(tile2lat(y + 1, z));\n var minx = lonToX(tile2lon(x, z));\n var maxx = lonToX(tile2lon(x + 1, z));\n var bb=[minx, miny, maxx, maxy]; \n \n function tile2lat(y, z) {\n var n = Math.PI - (2.0 * Math.PI * y) / Math.pow(2.0, z);\n return 180*(Math.atan(Math.sinh(n)))/Math.PI; //Lat in degrees\n }\n function tile2lon(x, z) {\n return (x / Math.pow(2.0, z) * 360.0 - 180); //Lon in degrees\n }\n function lonToX(lon) {\n return( lon * 20037508.34 / 180);\n }\n function latToY(lat) {\n return( 20037508.34*(Math.log(Math.tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180))/180);\n }\nreturn bb;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":787,"value":"md`### Step3:\nGet data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":40,"value":"slopeValues = readGeotiff2(wmsSlopeUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":229,"value":"slopeRange=findRange(slopeValues, \"excludezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":664,"value":"slopeDegValues = readGeotiff2(wmsSlopeDegUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":667,"value":"slopeDegRange = findRange(slopeDegValues, \"includezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":634,"value":"aspectValues = readGeotiff2(wmsAspectUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":672,"value":"aspectRange = findRange(aspectValues, \"includezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":687,"value":"aspectDegValues = readGeotiff2(wmsAspectDegUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":690,"value":"aspectDegRange = findRange(aspectDegValues, \"includezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":601,"value":"elevHillshadeValues = readGeotiff2(wmsElevationHillshadeUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":613,"value":"elevHillshadeRange= findRange(elevHillshadeValues, \"includezero\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":794,"value":"landCoverValues = readGeotiff2(wmsLandCoverUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":800,"value":"landCoverRange = findRange(landCoverValues, \"includezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":377,"value":"md`Compare elevation against a mapzen elevation tile to see if we got the right bounding box:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":356,"value":"wmsElevationUrlTif =\"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3ANone&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":774,"value":"elevationValues = readGeotiff2(wmsElevationUrlTif)","pinned":false,"hidden":false,"mode":"ojs"},{"id":776,"value":"elevationRange = findRange(elevationValues, \"includezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":410,"value":"mapzenValues = readGeotiff2(mapzenURL)","pinned":false,"hidden":false,"mode":"ojs"},{"id":406,"value":"mapzenRange= findRange(mapzenValues, \"includezero\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":505,"value":"//This is handy but return a cors error. One option is to add \"https://cors-anywhere.herokuapp.com/\" to the beginning of the url, but this limits the number of requests. Using the geotiff.fromArrayBuffer method instead (function readGeotiff2)\nfunction readGeotiff(url) {\n return geotiff.fromUrl(url)\n .then( tile => tile.getImage())\n .then( tileImage => tileImage.readRasters() )\n .then( tileReadRasters => tileReadRasters[0] )\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":614,"value":"function findRange(data, zeroRule) {\n var range=[9999999, -999999];\n for (var i=0; irange[1]){range[1]=data[i];}\n }\n return range;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":739,"value":"function readGeotiff2(url) {\n return fetch(url)\n .then(response => response.arrayBuffer())\n .then(arrayBuffer => geotiff.fromArrayBuffer(arrayBuffer))\n .then(tile => tile.getImage())\n .then(tileImage => tileImage.readRasters())\n .then(tileReadRasters => tileReadRasters[0]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":881,"value":"range = {\n if(map == \"3DEP slope pct png\"){ \n return [0,0];\n }\n if(map == \"3DEP slope pct geotiff\"){\n return slopeRange;\n }\n if(map == \"3DEP slope deg geotiff\"){\n return slopeDegRange\n }\n if(map == \"3DEP aspect deg geotiff\"){\n return aspectDegRange;\n }\n if(map == \"3DEP aspect-slope png\"){\n return [0,0];\n }\n if(map == \"3DEP aspect-slope geotiff\"){\n return aspectRange;\n }\n if(map == \"NLCD land cover png\"){\n return [0,0];\n }\n if(map == \"NLCD land cover geotiff\"){\n return landCoverRange;\n }\n if(map == \"3DEP elevation-hillshade png\"){\n return [0,0];\n }\n if(map == \"3DEP elevation-hillshade geotiff\"){\n return elevHillshadeRange;\n }\n if(map == \"3DEP elevation geotiff\"){\n return elevationRange;\n }\n if(map == \"Mapzen elevation geotiff\"){\n return mapzenRange;\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":123,"value":"md`### Step4: \nRender with Plotty (for geotiffs)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":46,"value":"function plottyRenderer(data, min, max, colorscale) {\n const canvas2 = DOM.canvas()\n var t0=performance.now();\n var plot = new plotty.plot({\n canvas: canvas2,\n data: data, width: 512, height: 512,\n domain: [min, max], colorScale: colorscale\n});\n //plot.setExpression('log');\n plot.render();\n var t1 =performance.now();\n var time = t1-t0 +\"ms\"\n console.log(\"Plotty render time= \"+time);\n return canvas2;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":112,"value":"md`### WMS urls:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":16,"value":"wmsSlopeUrlPng = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3ASlope%20Map&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Fpng&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":23,"value":"wmsSlopeUrlTif=\"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3ASlope%20Map&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":660,"value":"wmsSlopeDegUrlTif = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3ASlope%20Degrees&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":682,"value":"wmsAspectDegUrlTif = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3AAspect%20Degrees&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":553,"value":"wmsLandCoverUrlPng=\"https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2016_Land_Cover_L48/wms?service=WMS&request=GetMap&&Layers=NLCD%5F2016%5FLand%5FCover%5FL48&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&Format=image%2Fpng&width=512&height=512\"//\"https://www.mrlc.gov/geoserver/NLCD_Land_Cover/wms?SERVICE=WMS&REQUEST=GetMap&&Layers=NLCD%5FLand%5FCover&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&Format=image%2Fpng&width=512&height=512\"//https://www.mrlc.gov/data-services-page","pinned":false,"hidden":false,"mode":"ojs"},{"id":796,"value":"wmsLandCoverUrlTif=\"https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2016_Land_Cover_L48/wms?service=WMS&request=GetMap&&Layers=NLCD%5F2016%5FLand%5FCover%5FL48&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&Format=image%2Fgeotiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":547,"value":"wmsAspectUrlPng = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3AAspect%20Map&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Fpng&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":631,"value":"wmsAspectUrlTif = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3AAspect%20Map&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":590,"value":"wmsElevationHillshadeUrlPng = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3AHillshade%20Elevation%20Tinted&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Fpng&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":597,"value":"wmsElevationHillshadeUrlTif = \"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3AHillshade%20Elevation%20Tinted&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"//\"https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WMSServer?request=GetMap&service=WMS&version=1.1.1&layers=3DEPElevation%3AGreyHillshade_elevationFill&styles=default&srs=EPSG%3A3857&bbox=\"+bb[0]+\",\"+bb[1]+\",\"+bb[2]+\",\"+bb[3]+\"&format=image%2Ftiff&width=512&height=512\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":361,"value":"mapzenURL=\"https://s3.amazonaws.com/elevation-tiles-prod/geotiff/\"+tileCoords.z+\"/\"+tileCoords.x+\"/\"+tileCoords.y+\".tif\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":783,"value":"md`### Imports:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"d3 = require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"plotty = require('plotty@0.3.1/dist/plotty.js')","pinned":false,"hidden":false,"mode":"ojs"},{"id":26,"value":"geotiff = require(\"geotiff\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":278,"value":"import {select} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":828,"value":"import {number} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":429,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}