2b:null 2c:["$","title",null,{"children":"Slope and Aspect in Leaflet | Seth Lutske | Observable"}] 31:T649,html`

Slope and aspect are calculated from a digital elevation mode (DEM). The Mapbox Terrain RGB tilelayer is an RGB encoded DEM with global coverage. Each pixel's RGB value can be transformed into an elevation with the following function:

N:
Neighboring pixel to the north

W:
Neighboring pixel to the west

LatLng

E:
Neighboring pixel to the east

S:
Neighboring pixel to the south

`32:T5ba,md` \`\`\`javascript function height(R, G, B) { return -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1) } \`\`\` If we have the RGB value of the pixel at a \`LatLng\` coordinate, we can get its elevation. We can also get the elevation of the pixels that neighbor our coordinate in the four cardinal directions, as well as the distances between. We'll call distance between north and south \`dy\`, and the distance between east and west \`dx\`. Let's imagine a function not yet written called \`getElevation\`, which takes in a \`LatLng\` and outputs an elevation value in meters. \`\`\`javascript var point = map.project(LatLng, zoom) var N = map.unproject( { ...point, y: point.y - 1 }, zoom ), S = map.unproject( { ...point, y: point.y + 1 }, zoom ), E = map.unproject( { ...point, x: point.x + 1 }, zoom ), W = map.unproject( { ...point, x: point.x - 1 }, zoom ), eleN = getElevation(N), eleS = getElevation(S), eleE = getElevation(E), eleW = getElevation(W) var dx = map.distance(E, W), dy = map.distance(N, S) var dzdx = (eleE - eleW) / dx, dzdy = (eleN - eleS) / dy \`\`\` With those values, ***slope*** and ***aspect*** are calculated like this: \`\`\`javascript const slope = Math.atan( Math.sqrt( dzdx ** 2 + dzdy ** 2)) * (180 / Math.PI) const aspect = dx !== 0 // counterclockwise from east ? ( Math.atan2( dzdy, dzdx ) * (180 / Math.PI) + 180 ) % 360 : ( 90 * (dy > 0 ? 1 : -1) + 180 ) % 360 \`\`\` `33:T51d,md` ### Getting the Pixel ### Now that we have the tile image as a function of \`latlng\` and zoom, we need to be able to find the pixel within that image that is associated with that \`latlng\`: \`\`\`javascript var pixelCoords = { x: Math.floor(projectedPoint.x) - tileX * 256, y: Math.floor(projectedPoint.y) - tileY * 256 } \`\`\` When the image loads, we can get its content, paint it to a canvas, and then read the RGBA values of our pixel: \`\`\`javascript tileImage.onload = function(){ var c = document.createElement('canvas') c.width = c.height = 256 var demCtx = c.getContext('2d') demCtx.drawImage(tileImage, 0, 0) const imageData = demCtx.getImageData(0, 0, 256, 256) const pixelRGB = getRGBA(imageData, xyPositionOnTile.x, xyPositionOnTile.y) } \`\`\` Where \`pixelRGB\` is defined like this: \`\`\`javascript function getRGBA(imgData, x, y) { var index = y * imgData.width + x var i = index * 4 var d = imgData.data return { R: d[i], G: d[i+1], B: d[i+2], A: d[i+3] } } \`\`\` I chose this roundabout way of getting the RGB value because of the way I'm stashing my tile data once its initially loaded. More on that later. A simpler method would be: \`\`\`javascript const pixelRGB = demCtx.getImageData(xyPositionOnTile.x, xyPositionOnTile.y, 1, 1) \`\`\` `34:Ta23,html`

Latitude: Longitude: Zoom:

Tile coords:

Tile url: https://api.mapbox.com/v4/mapbox.terrain-rgb/7/67/35.pngraw?access_token=$\{Mapbox_Access_Token}

Tile previously saved:

Projected Point Coordinate (L.Map.project):

Pixel Position on Tile:

RGB Value at Pixel:

`35:T5bb,html`

Elevation:

Slope:

Aspect:

`36:T4c2,md`## Stashing the Data ## If you know you'll be needing slope and aspect data for many points in the same general area, stashing your already-fetched data may be a good idea. If you dig into the source code in the \`getTopography\` cell (found in the Appendix at the bottom), you'll notice that we don't fetch a tile image for every pixel request. The \`getRGBAfromProjectedPoint\` function saves the tile data to the \`tiles\` object for every newly fetched tile. If a \`LatLng\`'s tile is already saved, it will just read data from the \`tiles\` object. How you store the data is up to you. You can store each tile image as an \`ImageBitMap\`, which doesn't take up much space (about 40 *bytes*). But then you'll have to write the bitmap to a canvas and run \`getImageData\` every time you want a value, which might be slow. Alternatively, you can save the entire 256 x 256 x 4 \`Uint8ClampedArray\` returned by \`getImageData\`, which takes up almost 3 megabytes per tile. The initial load is a bit slow, but accessing the rgb values from there is very quick (this is what this notebook is currently doing). These are just two methods. If and how you store the tiles depends on your needs. -------- `37:T51d,md`## Conclusion ## Talking about accuracy when taking topographic measurements can get philosophical. Most topographical measurements are fractal in nature. [Like measuring the length of a coastline](https://en.wikipedia.org/wiki/Coastline_paradox), values depends on the size of your ruler. If its too big, we may miss mountains. If its too small, we may be measuring the slope of a snail's shell, not the hill it lives on. Do we need to know what direction an entire mountainside faces, or just one slope on the face of that mountainside? We need to choose the size of our ruler based on the needs of a given application. What grows out of this method is the ability to dynamically create DEM data for any map extent on the globe. Many applications that use gobs of DEM data require that it be downloaded ahead of time for good performance. This is outside the paradigm of internet slippymaps. But now it can be done. If you know the area that you need data for ahead of time, you can use [this formula](https://gist.github.com/mourner/8825883) or [xyz-affair](https://github.com/veltman/xyz-affair) to get all the tile coordinates of a certain zoom within a given \`LatLngBounds\`, and preload and save the tiles. This brings a new level of possibility to topographical analysis in leaflet. `38:T1f7a,PopupModifier = { // https://github.com/slutske22/leaflet-popup-modifier // Adding nametag labels to all popup-able leaflet layers const sourceTypes = ['Layer','Circle','CircleMarker','Marker','Polyline','Polygon','ImageOverlay','VideoOverlay','SVGOverlay','Rectangle','LayerGroup','FeatureGroup','GeoJSON'] sourceTypes.forEach( type => { L[type].include({ nametag: type.toLowerCase() }) }) // Adding new options to the default options of a popup L.Popup.mergeOptions({ removable: false, editable: false, }) // Modifying the popup mechanics L.Popup.include({ // modifying the _initLayout method to include edit and remove buttons, if those options are enabled // ---------------- Source code ---------------------------- // // original from https://github.com/Leaflet/Leaflet/blob/master/src/layer/Popup.js _initLayout: function () { var prefix = 'leaflet-popup', container = this._container = L.DomUtil.create('div', prefix + ' ' + (this.options.className || '') + ' leaflet-zoom-animated'); var wrapper = this._wrapper = L.DomUtil.create('div', prefix + '-content-wrapper', container); this._contentNode = L.DomUtil.create('div', prefix + '-content', wrapper); L.DomEvent.disableClickPropagation(wrapper); L.DomEvent.disableScrollPropagation(this._contentNode); L.DomEvent.on(wrapper, 'contextmenu', L.DomEvent.stopPropagation); this._tipContainer = L.DomUtil.create('div', prefix + '-tip-container', container); this._tip = L.DomUtil.create('div', prefix + '-tip', this._tipContainer); if (this.options.closeButton) { var closeButton = this._closeButton = L.DomUtil.create('a', prefix + '-close-button', container); closeButton.href = '#close'; closeButton.innerHTML = '×'; L.DomEvent.on(closeButton, 'click', this._onCloseButtonClick, this); } // ---------------- Source code ---------------------------- // // --------------- My additions --------------------------- // var nametag if (this.options.nametag){ nametag = this.options.nametag } else if (this._source) { nametag = this._source.nametag } else { nametag = "popup" } if (this.options.removable && !this.options.editable){ var userActionButtons = this._userActionButtons = L.DomUtil.create('div', prefix + '-useraction-buttons', wrapper); var removeButton = this._removeButton = L.DomUtil.create('a', prefix + '-remove-button', userActionButtons); removeButton.href = '#close'; removeButton.innerHTML = `Remove this ${nametag}`; this.options.minWidth = 110; L.DomEvent.on(removeButton, 'click', this._onRemoveButtonClick, this); } if (this.options.editable && !this.options.removable){ var userActionButtons = this._userActionButtons = L.DomUtil.create('div', prefix + '-useraction-buttons', wrapper); var editButton = this._editButton = L.DomUtil.create('a', prefix + '-edit-button', userActionButtons); editButton.href = '#edit'; editButton.innerHTML = 'Edit'; L.DomEvent.on(editButton, 'click', this._onEditButtonClick, this); } if (this.options.editable && this.options.removable){ var userActionButtons = this._userActionButtons = L.DomUtil.create('div', prefix + '-useraction-buttons', wrapper); var removeButton = this._removeButton = L.DomUtil.create('a', prefix + '-remove-button', userActionButtons); removeButton.href = '#close'; removeButton.innerHTML = `Remove this ${nametag}`; var editButton = this._editButton = L.DomUtil.create('a', prefix + '-edit-button', userActionButtons); editButton.href = '#edit'; editButton.innerHTML = 'Edit'; this.options.minWidth = 160; L.DomEvent.on(removeButton, 'click', this._onRemoveButtonClick, this); L.DomEvent.on(editButton, 'click', this._onEditButtonClick, this); } }, _onRemoveButtonClick: function (e) { this._source.remove(); L.DomEvent.stop(e); }, _onEditButtonClick: function (e) { //Needs to be defined first to capture width before changes are applied var inputFieldWidth = this._inputFieldWidth = this._container.offsetWidth - 2*19; this._contentNode.style.display = "none"; this._userActionButtons.style.display = "none"; var wrapper = this._wrapper; var editScreen = this._editScreen = L.DomUtil.create('div', 'leaflet-popup-edit-screen', wrapper) var inputField = this._inputField = L.DomUtil.create('div', 'leaflet-popup-input', editScreen); inputField.setAttribute("contenteditable", "true"); inputField.innerHTML = this.getContent() // ----------- Making the input field grow till max width ------- // inputField.style.width = inputFieldWidth + 'px'; var inputFieldDiv = L.DomUtil.get(this._inputField); // create invisible div to measure the text width in pixels var ruler = L.DomUtil.create('div', 'leaflet-popup-input-ruler', editScreen); let thisStandIn = this; // Padd event listener to the textinput to trigger a check this._inputField.addEventListener("keydown", function(){ // Check to see if the popup is already at its maxWidth // and that text doesnt take up whole field if (thisStandIn._container.offsetWidth < thisStandIn.options.maxWidth + 38 && thisStandIn._inputFieldWidth + 5 < inputFieldDiv.clientWidth){ ruler.innerHTML = inputField.innerHTML; if (ruler.offsetWidth + 20 > inputFieldDiv.clientWidth){ console.log('expand now'); inputField.style.width = thisStandIn._inputField.style.width = ruler.offsetWidth + 10 + 'px'; thisStandIn.update(); } } }, false) var inputActions = this._inputActions = L.DomUtil.create('div', 'leaflet-popup-input-actions', editScreen); var cancelButton = this._cancelButton = L.DomUtil.create('a', 'leaflet-popup-input-cancel', inputActions); cancelButton.href = '#cancel'; cancelButton.innerHTML = 'Cancel'; var saveButton = this._saveButton = L.DomUtil.create('a', 'leaflet-popup-input-save', inputActions); saveButton.href = "#save"; saveButton.innerHTML = 'Save'; L.DomEvent.on(cancelButton, 'click', this._onCancelButtonClick, this) L.DomEvent.on(saveButton, 'click', this._onSaveButtonClick, this) this.update(); L.DomEvent.stop(e); }, _onCancelButtonClick: function (e) { L.DomUtil.remove(this._editScreen); this._contentNode.style.display = "block"; this._userActionButtons.style.display = "flex"; this.update(); L.DomEvent.stop(e); }, _onSaveButtonClick: function (e) { var inputField = this._inputField; if (inputField.innerHTML.length > 0){ this.setContent(inputField.innerHTML) } else { alert('Enter something'); }; L.DomUtil.remove(this._editScreen); this._contentNode.style.display = "block"; this._userActionButtons.style.display = "flex"; this.update(); L.DomEvent.stop(e); // ---------------------End my additions --------------------------------------- // } }) return "Leaflet Removable and Editable Popup Plugin" }39:T46b,html` `3a:T11a9,mapWithColorpicker = { const mapoptions = { center: {lat: 61.1581999957237, lng: 7.959911974109774}, zoom: 15, maxZoom: 15, keyboard: false } const Mapbox_Access_Token = 'pk.eyJ1Ijoic2x1dHNrZTIyIiwiYSI6ImNqeGw1Y3BibDAybG4zeHFyaXl3OXVwZXUifQ.fZ_5Raq5z-DUpo2AK-bQHA' const map = L.map(document.getElementById('map0'),mapoptions) const RGBTerrain = colorPicker(`https://api.mapbox.com/v4/mapbox.terrain-rgb/{z}/{x}/{y}.pngraw?access_token=${Mapbox_Access_Token}`) RGBTerrain.addTo(map) function height(R, G, B) { return -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1) } function getElevation(latlng){ const rgb = RGBTerrain.getColor(latlng) return height(rgb[0], rgb[1], rgb[2]) } function getRGB(latlng){ const rgb = RGBTerrain.getColor(latlng) return { R: rgb[0], G: rgb[1], B: rgb[2] } } map.on('click', e => { // ---------- get dom elements ---------------------- const coordDiv = document.getElementById('coord') const Ndiv = document.getElementById('N') const Sdiv = document.getElementById('S') const Ediv = document.getElementById('E') const Wdiv = document.getElementById('W') // -------------------------------------------------- // ------- calculate pixels rgb and elevations ------ var point = map.project(e.latlng) var N = map.unproject( { ...point, y: point.y - 1 } ), S = map.unproject( { ...point, y: point.y + 1 } ), E = map.unproject( { ...point, x: point.x + 1 } ), W = map.unproject( { ...point, x: point.x - 1 } ), rgbN = getRGB(N), rgbS = getRGB(S), rgbE = getRGB(E), rgbW = getRGB(W), eleN = getElevation(N), eleS = getElevation(S), eleE = getElevation(E), eleW = getElevation(W) // -------------------------------------------------- // ----------- dxdz, dydz.slope and aspect ---------- var dx = map.distance(E, W), dy = map.distance(N, S) var dzdx = (eleE - eleW) / dx, dzdy = (eleN - eleS) / dy const slope = Math.atan( Math.sqrt( dzdx ** 2 + dzdy ** 2)) * (180 / Math.PI) const aspect = dx !== 0 ? 90 - Math.atan2( dzdy, dzdx ) * (180 / Math.PI) : 90 - 90 * (dy > 0 ? 1 : -1) // -------------------------------------------------- // ----------- popuplate and style DOM grid --------- Ndiv.style.backgroundColor = `rgb(${rgbN.R}, ${rgbN.G}, ${rgbN.B}` Sdiv.style.backgroundColor = `rgb(${rgbS.R}, ${rgbS.G}, ${rgbS.B}` Ediv.style.backgroundColor = `rgb(${rgbE.R}, ${rgbE.G}, ${rgbE.B}` Wdiv.style.backgroundColor = `rgb(${rgbW.R}, ${rgbW.G}, ${rgbW.B}` Ndiv.innerHTML = ` X: ${point.x.toFixed(0)}
Y: ${(point.y - 1).toFixed(0)}

R: ${rgbN.R}, G: ${rgbN.G}, B: ${rgbN.B})

Elevation: ${eleN.toFixed(1)} m ` Sdiv.innerHTML = ` X: ${point.x.toFixed(0)}
Y: ${(point.y + 1).toFixed(0)}

R: ${rgbS.R}, G: ${rgbS.G}, B: ${rgbS.B})

Elevation: ${eleS.toFixed(1)} m ` Ediv.innerHTML = ` X: ${(point.x + 1).toFixed(0)}
Y: ${point.y.toFixed(0)}

R: ${rgbE.R}, G: ${rgbE.G}, B: ${rgbE.B})

Elevation: ${eleE.toFixed(1)} m ` Wdiv.innerHTML = ` X: ${(point.x - 1).toFixed(0)}
Y: ${point.y.toFixed(0)}

R: ${rgbW.R}, G: ${rgbW.G}, B: ${rgbW.B})

Elevation: ${eleW.toFixed(1)} m ` coordDiv.innerHTML = ` Lat: ${e.latlng.lat.toFixed(3)}
Lng: ${e.latlng.lng.toFixed(3)}
dx: ${dx.toFixed(1)} m
dy: ${dy.toFixed(1)} m

Slope: ${slope.toFixed(1)}°
Aspect: ${aspect.toFixed(1)}° ` // -------------------------------------------------- const color = RGBTerrain.getColor(e.latlng) }) return map }3b:T6b2,mapWithDebug = { const mapoptions = { zoomControl: false, attributionControl: false, center: {lat: 61.1581999957237, lng: 7.959911974109774}, zoom: 7, maxZoom: 15, keyboard: false } const Mapbox_Access_Token = 'pk.eyJ1Ijoic2x1dHNrZTIyIiwiYSI6ImNqeGw1Y3BibDAybG4zeHFyaXl3OXVwZXUifQ.fZ_5Raq5z-DUpo2AK-bQHA' const map = L.map(document.getElementById('map'),mapoptions) const RGBTerrain = L.tileLayer(`https://api.mapbox.com/v4/mapbox.terrain-rgb/{z}/{x}/{y}.pngraw?access_token=${Mapbox_Access_Token}`) RGBTerrain.addTo(map) L.GridLayer.DebugCoords = L.GridLayer.extend({ createTile: function (coords) { var tile = document.createElement('div'); tile.innerHTML = `X: ${coords.x}, Y: ${coords.y}, Z: ${coords.z}` tile.classList.add('custom-tile') tile.style.outline = '1px solid white'; tile.style.color = 'white' tile.style.padding = '5px' tile.style.fontSize = '2em' return tile; } }); L.gridLayer.debugCoords = function(opts) { return new L.GridLayer.DebugCoords(opts); }; map.addLayer( L.gridLayer.debugCoords() ); const latInput = document.getElementById('latInput') const lngInput = document.getElementById('lngInput') const zoomInput = document.getElementById('zoomInput') map.on('click', e => { getTopography(e.latlng.lat, e.latlng.lng, map.getZoom(), "workspace") latInput.value = e.latlng.lat lngInput.value = e.latlng.lng zoomInput.value = map.getZoom() }) return map }3c:Ta17,esrimap = { const mapOptions = { center: {lat: 63.27850865468602, lng: -42.77447489000177}, zoom: 10, maxZoom: 18 } const map = L.map(document.getElementById('esrimap'),mapOptions) map.on('click', e => { }) let globalToken const authservice = 'https://www.arcgis.com/sharing/rest/oauth2/token' const client_id = 'r7Q4aCEr4j8qKbwN' const client_secret = 'a56716463cc742aba686ab9429ed87bb' async function getEsriToken(){ const url = `${authservice}?client_id=${client_id}&client_secret=${client_secret}&grant_type=client_credentials&expiration=${7200}` let access_token await fetch(url, { method: 'POST', }) .then(res => { return res.json() }) .then(res => { access_token = res.access_token globalToken = res.access_token }) .catch(err => console.error(err)) return access_token } async function createEsriLayer(layerType, options){ let token = await getEsriToken() var esriLayer = EsriLeaflet[layerType]({ ...options, token }) esriLayer.addTo(map) return esriLayer } var EsriSlopeDegreesMapLayer (async () => { EsriSlopeDegreesMapLayer = await createEsriLayer('imageMapLayer', { url: 'https://elevation.arcgis.com/arcgis/rest/services/WorldElevation/Terrain/ImageServer', opacity: 1, renderingRule: {rasterFunction: "Slope_Degrees_Map"} }) EsriSlopeDegreesMapLayer.addTo(map) })() function placeSlopeMarker(e){ var marker = L.marker(e.latlng) marker.bindPopup('Loading elevation data...', {autoPan: false, removable: true, minWidth: 400, maxWidth: 450}) marker.addTo(e.target) marker.openPopup() let timer = 0 let timerInterval = setInterval( () => {timer++},1) EsriLeaflet.identifyImage({ url: 'https://elevation.arcgis.com/arcgis/rest/services/WorldElevation/Terrain/ImageServer' }) .at(e.latlng) .token(globalToken) .setRenderingRule({rasterFunction: "Slope_Degrees"}) .run(function (error, identifyImageResponse, rawResponse) { if (error) { console.log(error); return; } clearInterval(timerInterval) var slope = identifyImageResponse.pixel.properties.value marker.getPopup().setContent(` Slope: ${slope}°
Time to get: ${timer}ms `) }) } map.on('click', placeSlopeMarker) return map }3d:T4066,getTopography = { const Mapbox_Access_Token = 'pk.eyJ1Ijoic2x1dHNrZTIyIiwiYSI6ImNqeGw1Y3BibDAybG4zeHFyaXl3OXVwZXUifQ.fZ_5Raq5z-DUpo2AK-bQHA' // ------------- Check if an image is in the cache (trustworthy?) ---------------------- function is_cached(src) { var image = new Image(); image.src = src; if (image.complete){ return 'True' } else { return 'False' } } // -------------------------------------------------------------------------------- var tiles = {} function loadImage(src) { return new Promise((resolve, reject) => { const img = new Image(); img.crossOrigin = "*" img.addEventListener("load", () => resolve(img)); img.addEventListener("error", err => reject(err)); img.src = src; }); }; async function getTopography(latArg, lngArg, zoomArg, filldom) { // ---------------- Add to Number.prototype ------------------------------------- if (typeof(Number.prototype.toRad) === "undefined") { Number.prototype.toRad = function() { return this * Math.PI / 180; } } // -------------------------------------------------------------------------------- // ------------ define dom nodes for reference ------------------------------------ const latInput = document.getElementById('latInput') const lngInput = document.getElementById('lngInput') const zoomInput = document.getElementById('zoomInput') const coordsOuput = document.getElementById('tilecoords') const urlOutput = document.getElementById('tileaddress') const tileImage = document.getElementById('tileImage') const cachedOutput = document.getElementById('cachedstatus') const projectedpointOutput = document.getElementById('projectedpoint') const pixelPositionOnTileOutput = document.getElementById('pointontile') const tilecanvasWorkspace = document.getElementById('tilecanvas') const rgbvalueatpixel = document.getElementById('rgbvalueatpixel') const pixelsample = document.getElementById('pixelsample') const coord15 = document.getElementById('coord15') const tileImageZoom15 = document.getElementById('tileImageZoom15') const specs15 = document.getElementById('specs15') const tileCanvasZoom15 = document.getElementById('tileCanvasZoom15') const N15 = document.getElementById('N15') const S15 = document.getElementById('S15') const E15 = document.getElementById('E15') const W15 = document.getElementById('W15') const finalslope = document.getElementById('finalslope') const finalaspect = document.getElementById('finalaspect') const finalelevation = document.getElementById('finalelevation') var imageOutput var tilecanvas if (filldom === 'workspace'){ imageOutput = tileImage tilecanvas = tilecanvasWorkspace } if (filldom === 'final'){ imageOutput = tileImageZoom15 tilecanvas = tileCanvasZoom15 } // -------------------------------------------------------------------------------- // ----------- take in arguments from map call, or grab input elements ------------- const lat = latArg || Number(latInput.value) const lng = lngArg || Number(lngInput.value) const zoom = zoomArg || Number(zoomInput.value) // -------------------------------------------------------------------------------- // --------------- get X Y Z coords of tile ---------------------------------------- function getTileCoords(lat, lng, zoom){ var projected = refmap.project({lat, lng}, zoom) const coords = { X: Math.floor(projected.x / 256), Y: Math.floor(projected.y / 256), Z: zoom } return coords } function getTileCoordsFromProjectedPoint(projected){ return { X: Math.floor(projected.x / 256), Y: Math.floor(projected.y / 256), Z: zoom } } // -------------------------------------------------------------------------------- // ------------ Get x y of point's pixel position within the tile ------------------ function pixelPositionOnTile(projectedPoint) { return { x: Math.floor(projectedPoint.x) - Math.floor(projectedPoint.x / 256) * 256, y: Math.floor(projectedPoint.y) - Math.floor(projectedPoint.y / 256) * 256 } } // -------------------------------------------------------------------------------- // -------- get rgba values from imageData for given xy position on tile --------- function getRGBA(imgData, x, y) { var index = y * imgData.width + x var i = index * 4 var d = imgData.data return { R: d[i], G: d[i+1], B: d[i+2], A: d[i+3] } } // -------------------------------------------------------------------------------- // ----- check if image was previously loaded and saved in tiles object ------------ function tileInStorage(coords) { if ( tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] ){ return 'True' } else { return 'False' } } // -------------------------------------------------------------------------------- // -------------- Calculate variables ----------------------------------------------- // const coords = getTileCoords(lat, lng, zoom) const projectedPoint = refmap.project({lat, lng}, zoom) const coords = getTileCoordsFromProjectedPoint(projectedPoint) const imageUrl = `https://api.mapbox.com/v4/mapbox.terrain-rgb/${coords.Z}/${coords.X}/${coords.Y}.pngraw?access_token=${Mapbox_Access_Token}` const formattedProjectPoint = {x: projectedPoint.x.toFixed(4), y: projectedPoint.y.toFixed(4)} const xyPositionOnTile = pixelPositionOnTile(projectedPoint) // -------------------------------------------------------------------------------- // -------------- Popuplate Information ----------------------------------------------- if (filldom === 'workspace'){ coordsOuput.innerHTML = `Tile coords: ${JSON.stringify(coords)}` urlOutput.innerHTML = `Tile url:

https://api.mapbox.com/v4/mapbox.terrain-rgb/${coords.Z}/${coords.X}/${coords.Y}.pngraw?access_token=$\{Mapbox_Access_Token}

` cachedOutput.innerHTML = `Tile previously saved: ${tileInStorage(coords)}` projectedpointOutput.innerHTML = `Projected Point Coordinate: (L.Map.project): ${JSON.stringify(formattedProjectPoint)}` pixelPositionOnTileOutput.innerHTML = `Pixel Position on Tile: ${JSON.stringify(xyPositionOnTile)}` } imageOutput.crossOrigin = "*" imageOutput.src = imageUrl // ----------------------------------------------------------------------------------- function onload2promise(obj){ return new Promise((resolve, reject) => { obj.onload = () => resolve(obj); obj.onerror = reject; }); } function saveTileData(){ var c = document.createElement('canvas') c.width = c.height = 256 var demCtx = c.getContext('2d') demCtx.drawImage(imageOutput, 0, 0) const imageData = demCtx.getImageData(0, 0, 256, 256) tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] = imageData } // ------------ Get RGBA from project point coordinates --------------------------- async function getRGBAfromProjectedPoint(projectedPoint){ const coords = getTileCoordsFromProjectedPoint(projectedPoint) const xy = pixelPositionOnTile(projectedPoint) const preloadedTileData = tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] // if image data already saved, just grab the rgba value from the pixel's position if ( preloadedTileData ){ return getRGBA(preloadedTileData, xy.x, xy.y) } // if not, get the image, get the rgba value, save the image as Uint16BitArray, // return rgba value in the meantime (will return rgba for every first time a tile is fetched) // every time after that, the image will have been stashed and we'll never get here: let rgba; await loadImage(imageUrl) .then( image => { var c = document.createElement('canvas') c.width = c.height = 256 var demCtx = c.getContext('2d') demCtx.drawImage(image, 0, 0) const imageData = demCtx.getImageData(0, 0, 256, 256) tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] = imageData rgba = getRGBA(imageData, xy.x, xy.y) // console.log('size of 1 tile:', sizeof(imageData.data), 'bytes') // createImageBitmap(image, 0, 0, 256, 256) // .then( ibm => { // console.log('bitmap size', sizeof(ibm), ibm ) // var testCanvas = document.createElement('canvas') // var testContext = testCanvas.getContext('2d') // testContext.drawImage(ibm, 0, 0) // const testData = testContext.getImageData(0, 0, 256, 256) // console.log('size of ImageData', sizeof(testData)) // }) }) return rgba } // -------------------------------------------------------------------------------- // ------------ Get elevation from RGBA--------------------- --------------------- function rgba2elevation(RGBA){ const { R, G, B } = RGBA return -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1) } // -------------------------------------------------------------------------------- // ----- Alternative: Fetch image, turn into blob, assign to image --------------- // if ( !tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] ){ // await fetch(imageUrl) // .then(response => response.blob()) // .then(image => { // const urlObject = URL.createObjectURL(image) // imageOutput.crossOrigin = "*" // imageOutput.src = urlObject // imageOutput.onload = function(){ // var c = document.createElement('canvas') // c.width = c.height = 256 // demCtx = c.getContext('2d') // demCtx.drawImage(imageOutput, 0, 0) // const imageData = demCtx.getImageData(0, 0, 256, 256) // tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] = imageData // } // // imageOutput.onload = saveTileData // // error here: why? // // await onload2promise(imageOutput, saveTileData) // }) // .catch(error => { // console.log(error) // // write function to fetch image with URL with different zoom // }) // } // ----------------------------------------------------------------------------------- // ----- Paint image data to canvas and get rgb values at pixelPositionOnTile ------- var demCtx; imageOutput.onload = async function(){ var c = document.createElement('canvas') c.width = c.height = 256 demCtx = c.getContext('2d') demCtx.drawImage(imageOutput, 0, 0) const imageData = demCtx.getImageData(0, 0, 256, 256) const pixelRGB = getRGBA(imageData, xyPositionOnTile.x, xyPositionOnTile.y) if (filldom === 'workspace'){ rgbvalueatpixel.innerHTML = `RGB Value at Pixel: ${JSON.stringify(pixelRGB)}` pixelsample.style.backgroundColor = `rgba(${pixelRGB.R}, ${pixelRGB.G}, ${pixelRGB.B}, ${pixelRGB.A})` pixelsample.innerHTML = `

R: ${pixelRGB.R}
G: ${pixelRGB.G}
B: ${pixelRGB.B}

` } if (filldom === 'final'){ // ------ CALCULATE SLOPE AND ASPECT AND SHOW IN DOM -------------------- var projectedN = { ...projectedPoint, y: projectedPoint.y - 4 }, projectedS = { ...projectedPoint, y: projectedPoint.y + 4 }, projectedE = { ...projectedPoint, x: projectedPoint.x + 4 }, projectedW = { ...projectedPoint, x: projectedPoint.x - 4 } var Nrgba = await getRGBAfromProjectedPoint(projectedN), Srgba = await getRGBAfromProjectedPoint(projectedS), Ergba = await getRGBAfromProjectedPoint(projectedE), Wrgba = await getRGBAfromProjectedPoint(projectedW), coordrgba = await getRGBAfromProjectedPoint(projectedPoint) var eleN = rgba2elevation(Nrgba), eleS = rgba2elevation(Srgba), eleE = rgba2elevation(Ergba), eleW = rgba2elevation(Wrgba), coordsElevation = rgba2elevation(coordrgba) var N = refmap.unproject( projectedN, 15 ), S = refmap.unproject( projectedS, 15 ), E = refmap.unproject( projectedE, 15 ), W = refmap.unproject( projectedW, 15 ) var dx = refmap.distance(E, W), dy = refmap.distance(N, S) var dzdx = (eleE - eleW) / dx, dzdy = (eleN - eleS) / dy const slope = Math.atan( Math.sqrt( dzdx ** 2 + dzdy ** 2)) * (180 / Math.PI) const aspect = dx !== 0 ? ( Math.atan2( dzdy, dzdx ) * (180 / Math.PI) + 180 ) % 360 : ( 90 * (dy > 0 ? 1 : -1) + 180 ) % 360 coord15.innerHTML = `Lat: ${lat.toFixed(4)}
Lng: ${lng.toFixed(4)}` specs15.innerHTML = ` Lat: ${lat.toFixed(4)}, Lng: ${lng.toFixed(4)}

Tile Coords: ${JSON.stringify(coords)}

R: ${pixelRGB.R}
G: ${pixelRGB.G}
B: ${pixelRGB.B}
` N15.style.backgroundColor = `rgb(${Nrgba.R},${Nrgba.G},${Nrgba.B})` S15.style.backgroundColor = `rgb(${Srgba.R},${Srgba.G},${Srgba.B})` E15.style.backgroundColor = `rgb(${Ergba.R},${Ergba.G},${Ergba.B})` W15.style.backgroundColor = `rgb(${Wrgba.R},${Wrgba.G},${Wrgba.B})` N15.innerHTML = ` R: ${Nrgba.R}
G: ${Nrgba.G}
B: ${Nrgba.B}
` S15.innerHTML = ` R: ${Srgba.R}
G: ${Srgba.G}
B: ${Srgba.B}
` E15.innerHTML = ` R: ${Ergba.R}
G: ${Ergba.G}
B: ${Ergba.B}
` W15.innerHTML = ` R: ${Wrgba.R}
G: ${Wrgba.G}
B: ${Wrgba.B}
` finalelevation.innerHTML = `Elevation: ${coordsElevation.toFixed(1)} meters` finalslope.innerHTML = `Slope: ${slope.toFixed(1)}°` finalaspect.innerHTML = `Aspect: ${aspect.toFixed(1)}°` // ------------------------------------------------------------------------------- } // save data to tiles object if ( !tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] ){ tiles[`X${coords.X}Y${coords.Y}Z${coords.Z}`] = imageData } } // ----------------------------------------------------------------------------------- // ------- Put a red dot on the canvas at the pixelPositionOnTile ------------------ var ctx = tilecanvas.getContext('2d') ctx.clearRect(0, 0, tilecanvas.width, tilecanvas.height); ctx.beginPath(); ctx.strokeStyle = "red"; ctx.fillStyle = "red"; ctx.arc(xyPositionOnTile.x, xyPositionOnTile.y, 5, 0, 2 * Math.PI); ctx.stroke(); ctx.fill() // ----------------------------------------------------------------------------------- } window.getTopography = getTopography return getTopography }3e:Tbb3,TopoLayer = { var mapboxAccessToken = 'pk.eyJ1Ijoic2x1dHNrZTIyIiwiYSI6ImNqeGw1Y3BibDAybG4zeHFyaXl3OXVwZXUifQ.fZ_5Raq5z-DUpo2AK-bQHA' var uniqueId = ( function(){ var lastId = 0 return function(){ return ++lastId } } )() var TopoLayer = L.GridLayer.extend({ // object to hold canvas contexts as they are created and updated _contexts: {}, // array to recieve worker objects when they get created _workers: [], // add worker initialization to beforeAdd Method beforeAdd: function (map) { map._addZoomLimit(this); // object to hold canvas contexts as they are created and updated this._contexts = {} // array to recieve worker objects when they get created this._workers = [] // grab topotype from options const { topotype } = this.options for (let i = 0; i < 16; i++){ var number = i < 9 ? `0${i+1}` : i + 1 this._workers[i] = new Worker(workers[topotype], {name: `Worker ${topotype} ${number}`}) this._workers[i].onmessage = e => this.updateTile(e, this) } }, // terminate all workers when layer is removed onRemove: function (map) { for (let i = 0; i < 16; i++){ this._workers[i].terminate() } this._removeAllTiles(); L.DomUtil.remove(this._container); map._removeZoomLimit(this); this._container = null; this._tileZoom = undefined; }, // createTile method required - creates a new tile of the gridlayer createTile: function(coords){ var tile = L.DomUtil.create('canvas', 'leaflet-tile') var size = this.getTileSize() tile.width = size.x tile.height = size.y var ctx = tile.getContext('2d') var demCtx var id = uniqueId() this._contexts[id] = ctx // define a new image element and its attributes var demImg = new Image() var { x, y, z } = coords demImg.crossOrigin = "*" demImg.src = `https://api.mapbox.com/v4/mapbox.terrain-rgb/${z}/${x}/${y}.pngraw?access_token=${mapboxAccessToken}` demImg.onload = function(){ var c = document.createElement('canvas') c.width = c.height = 256 demCtx = c.getContext('2d') demCtx.drawImage(demImg, 0, 0) redraw() } const redraw = () => { var data = { id: id } data.raster = demCtx.getImageData(0, 0, 256, 256) var workerIndex = (x + y) % this._workers.length this._workers[workerIndex].postMessage(data) } return tile }, updateTile: function(e, instance){ var ctx = instance._contexts[e.data.id] var imgData = ctx.createImageData(256, 256) var shades = e.data.shades imgData.data.set(shades) ctx.putImageData(imgData, 0, 0) } }) return TopoLayer }3f:T1a23,workers = { const worker4slopetext = new Blob([` self.dems = {} onmessage = function (e) { self.slopeshades = {} if (e.data === 'clear') { self.slopeshades = {} return } if (e.data.raster) { const { data } = e.data.raster self.slopeshades[e.data.id] = raster2slopes(data) self.shades = shading(self.slopeshades[e.data.id]) } postMessage({ id: e.data.id, message: 'from worker', ele: self.slopeshades[e.data.id], shades: self.shades }) } function raster2dem (data){ const dem = new Int16Array(256 * 256) var x, y, i, j function height (R, G, B) { return -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1) } for (x = 0; x < 256; x++) { for (y = 0; y < 256; y++){ i = x + y * 256 j = i * 4 dem[i] = height( data[j], data[j + 1], data[j + 2] ) } } return dem } function raster2slopes(raster){ const dem = raster2dem(raster) const slopes = new Float32Array( 256 * 256 ) var x, y, dx, dy, i, j for (x = 1; x < 255; x++){ for (y = 1; y < 255; y++){ i = y * 256 + x dx = ((dem[i - 255] + 2 * dem[i + 1] + dem[i + 257]) - (dem[i - 257] + 2 * dem[i - 1] + dem[i + 255])) / 8; dy = ((dem[i + 255] + 2 * dem[i + 256] + dem[i + 257]) - (dem[i - 257] + 2 * dem[i - 256] + dem[i - 255])) / 8; slopes[i] = Math.atan( Math.sqrt( dx * dx + dy * dy ) ) } } /** Shameless Hack: * When calculating slope, we can't get values * that are on the edge of a tile, as we can't * get their neighbors. Their neighbors are on * a different tile. Rather than trying to coordinate * data between tiles, this loop takes the pixel * 2 pixels deep from each edge and copies it to * the edge pixel. The hack is 1px wide and barely * visible */ for (x = 0; x < 256; x++){ for (y = 0; y < 256; y++){ i = y * 256 + x if (x === 0){ j = y * 256 + x + 1 slopes[i] = slopes[j] } if (x === 255){ j = y * 256 + x - 1 slopes[i] = slopes[j] } if (y === 0){ j = (y + 1) * 256 + x slopes[i] = slopes[j] } if (y === 255){ j = (y - 1) * 256 + x slopes[i] = slopes[j] } } } return slopes } function shading(slopes){ var px = new Uint8ClampedArray( 256 * 256 * 4 ) for (let i = 0; i < slopes.length; i++){ var color = slopes[i] * 255 / (Math.PI / 2) px[4*i + 0] = color px[4*i + 1] = color px[4*i + 2] = color px[4*i + 3] = 255 } return px } `], {type: "text/javascript"}) const worker4slope = URL.createObjectURL(worker4slopetext); const worker4aspecttext = new Blob([` self.slopeaspects = {} onmessage = function (e) { if (e.data.raster) { const { data } = e.data.raster self.slopeaspects[e.data.id] = raster2slopeaspect(data) self.shades = shading(self.slopeaspects[e.data.id].slopes, self.slopeaspects[e.data.id].aspects) } postMessage({ id: e.data.id, message: 'from worker', ele: self.slopeaspects[e.data.id], shades: self.shades }) } function raster2dem(data){ const dem = new Int16Array(256 * 256) var x, y, i, j function height (R, G, B) { return -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1) } for (x = 0; x < 256; x++) { for (y = 0; y < 256; y++){ i = x + y * 256 j = i * 4 dem[i] = height( data[j], data[j + 1], data[j + 2] ) } } return dem } function raster2slopeaspect(raster){ const dem = raster2dem(raster) const aspects = new Float32Array( 256 * 256 ) const slopes = new Float32Array( 256 * 256 ) var x, y, dx, dy, i for (x = 1; x < 255; x++){ for (y = 1; y < 255; y++){ i = y * 256 + x dx = ((dem[i - 255] + 2 * dem[i + 1] + dem[i + 257]) - (dem[i - 257] + 2 * dem[i - 1] + dem[i + 255])) / 8; dy = ((dem[i + 255] + 2 * dem[i + 256] + dem[i + 257]) - (dem[i - 257] + 2 * dem[i - 256] + dem[i - 255])) / 8; aspects[i] = dx !== 0 ? 90 - Math.atan2( dy, -dx ) * (180 / Math.PI) : 90 - 90 * (dy > 0 ? 1 : -1) slopes[i] = Math.atan( Math.sqrt( dx * dx + dy * dy ) ) * 180 / Math.PI } } return {slopes, aspects} } function shading(slopes, aspects){ var px = new Uint8ClampedArray( 256 * 256 * 4 ) for (let i = 0; i < aspects.length; i++){ var rgb = hypsotint(slopes[i], aspects[i]) px[4*i + 0] = rgb.R px[4*i + 1] = rgb.G px[4*i + 2] = rgb.B px[4*i + 3] = 255 } return px } var colorsrgb = [ { R: 154, G: 251, B: 012 }, { R: 000, G: 173, B: 067 }, { R: 000, G: 104, B: 192 }, { R: 108, G: 000, B: 163 }, { R: 202, G: 000, B: 156 }, { R: 255, G: 085, B: 104 }, { R: 255, G: 171, B: 071 }, { R: 244, G: 250, B: 000 }, { R: 154, G: 251, B: 012 }, ] var breakpoints = [ 0, 22.5, 67.5, 112.5, 157.5, 202.5, 247.5, 292.5, 337.5, 360] function hypsotint(slope, aspect){ let correctedAspect = aspect < 0 ? 360 + aspect % 360 : aspect > 360 ? aspect % 360 : aspect for (let i = 0; i < breakpoints.length - 1; i++){ if (breakpoints[i] < correctedAspect && correctedAspect <= breakpoints[i + 1]){ if (slope > 1){ return colorsrgb[i] } return { R: 120, G: 120, B: 120 } } } return { R: 154, G: 251, B: 012 } } `], {type: "text/javascript"}) const worker4aspect = URL.createObjectURL(worker4aspecttext); const workers = { slope: worker4slope, aspect: worker4aspect } return workers }40:T796,sizeof = { /* sizeof.js A function to calculate the approximate memory usage of objects Created by Kate Morley - http://code.iamkate.com/ - and released under the terms of the CC0 1.0 Universal legal code: http://creativecommons.org/publicdomain/zero/1.0/legalcode */ /* Returns the approximate memory usage, in bytes, of the specified object. The * parameter is: * * object - the object whose size should be determined */ function sizeof(object){ // initialise the list of objects and size var objects = [object]; var size = 0; // loop over the objects for (var index = 0; index < objects.length; index ++){ // determine the type of the object switch (typeof objects[index]){ // the object is a boolean case 'boolean': size += 4; break; // the object is a number case 'number': size += 8; break; // the object is a string case 'string': size += 2 * objects[index].length; break; // the object is a generic object case 'object': // if the object is not an array, add the sizes of the keys if (Object.prototype.toString.call(objects[index]) != '[object Array]'){ for (var key in objects[index]) size += 2 * key.length; } // loop over the keys for (var key in objects[index]){ // determine whether the value has already been processed var processed = false; for (var search = 0; search < objects.length; search ++){ if (objects[search] === objects[index][key]){ processed = true; break; } } // queue the value to be processed if appropriate if (!processed) objects.push(objects[index][key]); } } } // return the calculated size return size; } return sizeof }41:[] 2f:{"id":"3b063f002ab455ab","slug":"slope-as-a-function-of-latlng-in-leaflet","title":"Slope and Aspect in Leaflet","thumbnail":"fe88b5d3008d5b0872d6cd5c6bfb818f240cab86396ae01d6b48fa9ba9534870","access_level":"public","create_time":"$D2020-05-08T15:29:55.064Z","update_time":"$D2020-10-19T17:55:31.674Z","view_time":null,"first_public_time":"$D2020-05-09T03:09:34.918Z","fork":false,"stars":9,"starred":false,"owner":{"type":"team","id":"e2585cdca261ab5e","login":"slutske22","name":"Seth Lutske","avatar_url":"https://avatars-next.observableusercontent.com/940b7490fa472faafb71f8cdaac26c287d7d2d7ff18c0501e01788756090c541","tier":"starter_2024"},"version":1470,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"4461d9138c16adbb","login":"slutske22","name":"Seth Lutske","avatar_url":"https://avatars-next.observableusercontent.com/940b7490fa472faafb71f8cdaac26c287d7d2d7ff18c0501e01788756090c541","tier":null}],"members":[{"type":"individual","id":"4461d9138c16adbb","login":"slutske22","name":"Seth Lutske","avatar_url":"https://avatars-next.observableusercontent.com/940b7490fa472faafb71f8cdaac26c287d7d2d7ff18c0501e01788756090c541","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":3,"collection":null,"type":"notebook","body":{"title":"Slope and Aspect in Leaflet","cells":[{"id":0,"value":"md`# Slope and Aspect in Leaflet`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":638,"value":"html`\n\n\n

\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":35,"value":"md`### ***Part 1: Slope and Aspect as a function of \\`LatLng\\`*** ###\n*This notebook demonstrates how to use Mapbox RGB Terrain tiles to get slope or aspect as a function of \\`LatLng\\`. Every Mapbox RGB tile is packed with 256 x 256 data points. We can harvest data from a \\`LatLng\\`'s associated hi-resolution tile, and we can save tile data in a user session to create a reliable and quick data source for the whole globe.*\n\n*Skip to **[Part 2: Visualizing Topography](https://observablehq.com/@slutske22/slope-and-aspect-in-leaflet)** for less math and more colors!*\n\n#### ⏰ 12 minute read ####\n---------\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":707,"value":"md`## Introduction ##`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":685,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":676,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":731,"value":"md`### getElevation(\\`LatLng\\`): The Quick and Dirty Way ###\nThe [tilelayer-colorpicker](https://github.com/frogcat/leaflet-tilelayer-colorpicker) plugin is a quick way to get the rgb value of a pixel at a certain \\`LatLng\\`. So we can quickly get our elevation like this:\n\n\\`\\`\\`javascript\nconst terrainRGB = L.tileLayer.colorPicker(Mapbox_Terrain_RGB_URL).addTo(map)\n\nfunction getElevation(latlng){\n const rgb = terrainRGB.getColor(latlng)\n return height(rgb[0], rgb[1], rgb[2])\n}\n\\`\\`\\`\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":808,"value":"md`\n### Test it Out ###\nLet's set up a map that will run these calculations for us. We can see the results in a live version of our pixel grid. Click around the map:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":751,"value":"html`\n

\n\n

\n\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":817,"value":"md`### What's wrong here? ###\nYou'll quickly realize a glaring issue. Our results change drastically based on our zoom value. The lower the zoom level, the farther apart the neighboring pixels. For example, at a zoom level of 4, neighboring pixels can be up tp 20km apart! We'd be calculating the slope and aspect from points 10km in each direction. We can't trust this method at lower zoom levels. But we can choose better data.\n\n-------\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":821,"value":"md` ## Handpicking Tiles ##\n\nTo get better accuracy, we need to get the pixel values of a \\`LatLng\\` and its neighbors from tiles at higher zoom levels. Mapbox RGB Terrain has novel tiles for zoom levels of up to 15. At zoom level 15, the distance between pixels is on the order of ~5 meters. This is a much more reliable distance to use to do our calculations. But for any given \\`LatLng\\`, we need to get the tile that corresponds to that coordinate, at a zoom level of 15, regardless of the current zoom level of the map. \n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":55,"value":"md`\n### Getting the Tile ###\nTo get the tile associated with a certain \\`LatLng\\` and zoom level, we can use the following formula:\n\n\\`\\`\\`javascript\nconst projectedPoint = refmap.project({lat, lng}, zoom)\n\nfunction getTileCoord(projectedPoint, zoom) {\n return {\n X: Math.floor(projectedPoint.x / 256),\n Y: Math.floor(projectedPoint.y / 256),\n Z: zoom\n }\n}\n\\`\\`\\`\n\nThen we can use our coordinates to get the image URL of the tile, and create an \\`\\` element to load it:\n\n\\`\\`\\`javascript\nconst { X, Y, Z } = getTileCoord(lat, lng, zoom)\nconst tileURL = \n \\`https://api.mapbox.com/v4/mapbox.terrain-rgb/\\${X}/\\${Y}/\\${Z}.pngraw?access_token=\\${Mapbox_Access_Token}\\`\n\nconst tileImage = new Image()\ntileImage.crossOrigin = \"*\"\ntileImage.src = tileURL\n\\`\\`\\`\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":376,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":319,"value":"md`\n### Test it Out ###\nLet's create a workspace where we can test this out. Using either the inputs or the map, the following cells will be able to tell you the coordinates and url of the a slippymap tile, as well as the \\`xy\\` value of the \\`LatLng\\` we're asking for. We can also get its RGB value. We'll also keep track of whether not a tile has already been accessed and saved.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":59,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":310,"value":"md`\n*** For convenience, here's a map that runs the function when clicked: ***`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":221,"value":"html`\n\n\n

\n\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1281,"value":"md`\nNow we just need to get our \\`N\\`, \\`S\\`, \\`W\\`, and \\`W\\` pixels, and run them through the same process that we just described to get their RGB values.\n\n-----------\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":902,"value":"md`## Balancing Quality and Coverage ##\n\nWe can execute this whole process while explicitly choosing high zoom level tiles for high resolution data. Although Mapbox RGB Terrain goes up to zoom level 15, I've found that if we back off to zoom levels of 14 or 13, we can still have good accuracy, but we cover more ground with each tile. We'll stick with zoom 15 for this next example. Click on the map below and you'll see that regardless of what zoom level the map is at, we'll be fetching the RGB tile with a zoom level of 15. We don't even need to add the RGB layer to the map, so we'll choose something more interesting:\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":914,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1277,"value":"md`\n---------\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1034,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1107,"value":"md`## Bonus Shortcut: Use Esri ###\n\nThis discussion would not be complete without mention of [esri-leaflet](https://esri.github.io/esri-leaflet/). We can configure esri-leaflet's \\`query\\` function on their [terrain](https://www.arcgis.com/home/item.html?id=58a541efc59545e6b7137f961d7de883) dataset to give slope or aspect as a function of \\`\\`LatLng\\`\\`:\n\n\\`\\`\\`javascript\nL.esri.identifyImage({\n url: 'https://elevation.arcgis.com/arcgis/rest/services/WorldElevation/Terrain/ImageServer'\n})\n .token(token)\n .at(LatLng)\n .setRenderingRule({rasterFunction: \"Slope_Degrees\"}) // or {rasterFunction: \"Aspect\"}\n .run(function (error, identifyImageResponse, rawResponse) {\n\n if (error) { console.log(error); return; }\n\n var slope = identifyImageResponse.pixel.properties.value\n\n })\n\\`\\`\\`\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1241,"value":"html`\n\n\n

Click the map to run this function:

\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1270,"value":"md`\nEsri calculates slope and aspect server-side on the fly from their DEM, and sends back the value. Sometimes their response is quick, sometimes its not. If you only need topographic data for a few points, this is probably OK. But in roughly the same time it takes to ask Esri for one datapoint for one \\`LatLng\\`, we can download and calculate elevation, slope, and aspect values for 254 * 254 points with a Mapbox RGB tile.\n\n---------\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1210,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1306,"value":"md`*Head on over to **[Part 2: Visualizing Topography](https://observablehq.com/@slutske22/slope-and-aspect-in-leaflet)** for less math and more visuals.*`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1332,"value":"html`\n\n

Seth Lutske

San Diego, California

\n\n${fab`github`}\n${fab`linkedin`}\n${fab`facebook`}\n${fa`envelope`}\n\n

\n\n\n `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":12,"value":"md`

\n\n--- \n\n## Appendix\n\n*Cells used to implement this notebook lie below. Click the left margin to toggle them open.*`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":15,"value":"L = require('leaflet')\n\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":22,"value":"html``","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1250,"value":"EsriLeaflet = require('esri-leaflet')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1256,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1262,"value":"$39","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":773,"value":"colorPicker = {\n L.TileLayer.ColorPicker = L.TileLayer.extend({\n options: {\n crossOrigin: \"anonymous\"\n },\n getColor: function(latlng) {\n var size = this.getTileSize();\n var point = this._map.project(latlng, this._tileZoom).floor();\n var coords = point.unscaleBy(size).floor();\n var offset = point.subtract(coords.scaleBy(size));\n coords.z = this._tileZoom;\n var tile = this._tiles[this._tileCoordsToKey(coords)];\n if (!tile || !tile.loaded) return null;\n try {\n var canvas = document.createElement(\"canvas\");\n canvas.width = 1;\n canvas.height = 1;\n var context = canvas.getContext('2d');\n context.drawImage(tile.el, -offset.x, -offset.y, size.x, size.y);\n return context.getImageData(0, 0, 1, 1).data;\n } catch (e) {\n return null;\n }\n }\n });\n L.tileLayer.colorPicker = function(url, options) {\n return new L.TileLayer.ColorPicker(url, options);\n };\n \n return L.tileLayer.colorPicker\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1397,"value":"md`--------`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":641,"value":"bannermap = {\n\n const mapOptions = {\n zoomControl: false,\n attributionControl: false,\n center: {lat: 21.189293538320758, lng: -157.1152210235596},\n zoom: 15,\n maxZoom: 15\n }\n\n const Mapbox_Access_Token = 'pk.eyJ1Ijoic2x1dHNrZTIyIiwiYSI6ImNqeGw1Y3BibDAybG4zeHFyaXl3OXVwZXUifQ.fZ_5Raq5z-DUpo2AK-bQHA'\n\n const map = L.map(document.getElementById('map3'),mapOptions)\n\n const RGBTerrain = L.tileLayer(`https://api.mapbox.com/v4/mapbox.terrain-rgb/{z}/{x}/{y}.pngraw?access_token=${Mapbox_Access_Token}`)\n RGBTerrain.addTo(map)\n\n map.on('click', () => {console.log(map.getCenter(), map.getZoom())} )\n\n const aspectlayer = new TopoLayer({ topotype: \"slope\" })\n aspectlayer.addTo(map)\n\n return map\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":767,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":339,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":344,"value":"refmap = {\n \n const mapoptions = {\n zoomControl: false,\n attributionControl: false,\n center: {lat: 61.1581999957237, lng: 7.959911974109774},\n zoom: 7,\n maxZoom: 15\n }\n\n const map = L.map(document.getElementById('refmap'), mapoptions)\n return map\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":944,"value":"mapZoom15 = {\n\n\n const mapOptions = {\n center: {lat: 20.77092112928836, lng: -156.32285476887617},\n zoom: 10,\n maxZoom: 18\n }\n\n const map = L.map(document.getElementById('map15'),mapOptions)\n \n const Mapbox_Access_Token = 'pk.eyJ1Ijoic2x1dHNrZTIyIiwiYSI6ImNqeGw1Y3BibDAybG4zeHFyaXl3OXVwZXUifQ.fZ_5Raq5z-DUpo2AK-bQHA'\n \n // map.on('click', e => console.log(e.latlng, map.getZoom()))\n \n map.on('click', e => {\n \n getTopography(e.latlng.lat, e.latlng.lng, 15, 'final')\n \n })\n \nvar Stamen_Terrain = L.tileLayer('https://stamen-tiles-{s}.a.ssl.fastly.net/terrain/{z}/{x}/{y}{r}.{ext}', {\n\tattribution: 'Map tiles by Stamen Design, CC BY 3.0 — Map data © OpenStreetMap contributors',\n\tsubdomains: 'abcd',\n\tminZoom: 0,\n\tmaxZoom: 18,\n\text: 'png'\n}).addTo(map)\n \n return map\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1244,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1400,"value":"md`----------`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":299,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":662,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":665,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1447,"value":"md`----------`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":165,"value":"{\n \n Number.prototype.toRad = function() {\n return this * Math.PI / 180;\n }\n \n return \"Adding Number.prototype.toRad\"\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1122,"value":"$40","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1349,"value":"import { style as faStyle, fa } from \"@airbornemint/fontawesome\"\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1357,"value":"faStyle({regular: true, solid: true, brands: true})\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1375,"value":"function faCommon(kind) {\n return function(s) {\n return ``;\n };\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1378,"value":"far = faCommon(\"r\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1372,"value":"fas = faCommon(\"s\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1381,"value":"fab = faCommon(\"b\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W41","stdlib":"1","titling":"auto"}} 30:"$2f:body"