2d:null 2e:["$","title",null,{"children":"Northern Colorado Land Use Changes | Graham Stewart | Observable"}] 33:T6c9,

Northern Colorado Land Use Changes

2017 to 2021

The map below allows you to toggle between 2017 and 2021 and view the differences in land use in that 4 year period. All land use data comes from the 10m Annual Land Use Survey which is created by feeding ESA Sentinal-2 imagery into a Deep Neural Network run on Microsoft's Planetary Computer. Each 10x10m square of the grid is classified into a 9-category system and these are displayed color-coded on the map.

As you scroll the map, the Sankey chart below it will update to reflect how the different land categories have changed.

Known issues

The selected dataset only includes the MGRS 13T zone - as you scroll south it'll stop finding map tiles. I think I could override the getTileURL function to support this, but it's beyond the scope of this project.
The Sankey calculation is done clientside on the fetched land use tiles, this makes it easy to implenent as live, but it can be a little slow on some computers and also means that it doesn't have access to more granular zoom levels and small changes could be missed
Computation of acreage changes is done by calculating the total area of the Map Boundary and then assuming each of the pixels is equal area. The nature of the EPSG:3857 Web Mercator projection means that this assumption becomes increasingly inaccurate as more degrees of latitude are included. For the scales shown here, this is usable, but a better calculation would be necessary to calculate ice-loss in Greenland.

34:T74a, viewof map = { // Basic Leaflet / Observable code from https://observablehq.com/@coreydockser/learning-leaflet console.log("initializing map"); let parent = DOM.element('div', { style: `width:${width}px;height:${width/2}px` }); // Make sure Observable renders the parent container so Leaflet can figure out it's size yield parent; // Create a map centered on Northeastern's campus let myMap = L.map(parent, {drawControl: true}).setView([ 40.163657614276126, -105.20164489746095], 11); // higher numbers are more zoomed in // this sets the map into the value of this cell so we can access it later parent.value = myMap; // This function let mapChanged = function(ev) { console.log('fire'); parent.dispatchEvent(new CustomEvent('input')); }; // myMap.addEventListener('load', mapChanged); // Promises.delay(1000, mapChanged()); myMap.addEventListener('moveend', mapChanged); let osmTileLayer = L.tileLayer( 'http://{s}.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png', { attribution: '© OpenStreetMap' }).addTo(myMap ); osmTileLayer.name="base"; let pcLayer = L.tileLayer( 'https://planetarycomputer.microsoft.com/api/data/v1/item/tiles/{z}/{x}/{y}.png?scale=1&collection=io-lulc-9-class&item=13T-2021&assets=data&colormap_name=io-lulc-9-class', { attribution: 'Microsoft Planetary Computer' }).setOpacity(0.3).addTo( myMap ); pcLayer.name="2021"; let pcLayer2 = L.tileLayer( 'https://planetarycomputer.microsoft.com/api/data/v1/item/tiles/{z}/{x}/{y}.png?scale=1&collection=io-lulc-9-class&item=13T-2017&assets=data&colormap_name=io-lulc-9-class', { attribution: 'Microsoft Planetary Computer' }).setOpacity(0.3).addTo( myMap ); pcLayer.setOpacity(0); pcLayer2.name="2017"; setTimeout(mapChanged, 50); }35:T14d1,async function performAnalysis() { let result={}; console.log(map.getCenter()); console.log(map.getZoom()); console.log("starting analyze"); const NO_DATA =0; //255,255,255 const WATER=1; // 65,155,223 const TREES=2; // 57,125, 73 const FLOODED_VEGETATION=3; //122,135,198 const CROPS=4; //228,150, 53 const BUILT_AREA = 5; //196, 40, 27 const BARE_GROUND=6; //165,155,143 const SNOW_ICE=7; //168,235,255 const CLOUDS=8; // 97, 97, 97 const RANGELAND=9; //227,226,195 const FRIENDLY_NAMES = ["No Data","Water","Trees","Flooded Vegetation","Crops","Built Area","Bare Ground","Snow/Ice","Clouds","Rangeland"]; const SQM_PER_ACRE = 4046.86; //ugh, why no metric in this country const WGS_EARTH_RADIUS = 6378137.0; // used to calculate areas // ripped from the Leaflet-draw library because it wasn't happy being included here and is total ovekill let geodesicArea = function (latLngs) { var pointsCount = latLngs.length, area = 0.0, d2r = Math.PI / 180, p1, p2; if (pointsCount > 2) { for (var i = 0; i < pointsCount; i++) { p1 = latLngs[i]; p2 = latLngs[(i + 1) % pointsCount]; area += ((p2.lng - p1.lng) * d2r) * (2 + Math.sin(p1.lat * d2r) + Math.sin(p2.lat * d2r)); } area = area * (WGS_EARTH_RADIUS*WGS_EARTH_RADIUS) / 2.0; } return Math.abs(area); } // We can determine the use type entirely by the red pixels as each item in the symbology varies function getUseType(redPixel) { switch (redPixel) { case 65: return WATER; case 57: return TREES; case 122: return FLOODED_VEGETATION; case 228: return CROPS; case 196: return BUILT_AREA; case 165: return BARE_GROUND; case 168: return SNOW_ICE; case 97: return CLOUDS; case 227: return RANGELAND; } return redPixel; // console.log("Undefined r of "+ redPixel); } // This is a custom forked version of the LeafletImage library that will let me render a single layer // instead of the map as it shows up on screen. let complete=false; let ls2017only = function(l){ return l.name=="2017";} LeafletImage(map, ls2017only, function(err,canvas2017) { let ls2021Only = function(l){ return l.name=="2021";} LeafletImage(map, ls2021Only, function(err,canvas2021) { var dimensions = map.getSize(); const image2017 = canvas2017.getContext("2d").getImageData(0,0,dimensions.x,dimensions.y); const image2021 = canvas2021.getContext("2d").getImageData(0,0,dimensions.x,dimensions.y); result.totalMapAreaSqM = geodesicArea( [map.getBounds().getNorthWest(), map.getBounds().getNorthEast(), map.getBounds().getSouthEast(), map.getBounds().getSouthWest() ]); result.acresPerPixel = result.totalMapAreaSqM/ SQM_PER_ACRE / (dimensions.x * dimensions.y); result.pixelsChangedDict=[]; result.totalPixels=0; result.pixelsChanged=0; for (let i = 0; i < image2017.data.length; i += 4) { result.totalPixels++; let pixelUse2017 =getUseType(image2017.data[i]); let pixelUse2021 =getUseType(image2021.data[i]); if (pixelUse2017!=pixelUse2021) { result.pixelsChanged++; } if (!result.pixelsChangedDict[pixelUse2017 + "|" + pixelUse2021]) { result.pixelsChangedDict[pixelUse2017 + "|" + pixelUse2021]=0; } result.pixelsChangedDict[pixelUse2017 + "|" + pixelUse2021]++; } result.sankey=[]; for (let src=0;src1000) { result.sankey.push({source:FRIENDLY_NAMES[src] , target: FRIENDLY_NAMES[dest]+" ", value:(result.pixelsChangedDict[src + "|" + dest]*result.acresPerPixel) }); } } } console.log("done analyzing"); complete=true; console.log(result); }); }); // THis is horribly hacky, but stops us returning until everything is loaded while (!complete) { await new Promise(r => setTimeout(r, 100)); } return result; }36:[] 31:{"id":"1843ad814e3fe3bc","slug":"colorado-land-use-changes","title":"Northern Colorado Land Use Changes","thumbnail":"9f10e59e14d93eeb8bc34a39fc178840ab8ac6385067457763c50815cbca8a64","access_level":"public","create_time":"$D2022-10-08T03:40:25.165Z","update_time":"$D2022-11-13T17:17:11.685Z","view_time":null,"first_public_time":"$D2022-11-13T16:26:21.283Z","has_unpublished_edits":null,"fork":false,"stars":0,"starred":false,"owner":{"type":"team","id":"1b919c6c0c8a2266","login":"grahamsz","name":"Graham Stewart","avatar_url":"https://avatars-next.observableusercontent.com/cbf473b0877318cca19e6fdd40b797886f76e8defaa989ca59a07afb535b0a3f","tier":"starter_2024"},"version":668,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"713490e721b1ecbd","login":"grahamsz","name":"Graham Stewart","avatar_url":"https://avatars-next.observableusercontent.com/cbf473b0877318cca19e6fdd40b797886f76e8defaa989ca59a07afb535b0a3f","tier":null}],"members":[{"type":"individual","id":"713490e721b1ecbd","login":"grahamsz","name":"Graham Stewart","avatar_url":"https://avatars-next.observableusercontent.com/cbf473b0877318cca19e6fdd40b797886f76e8defaa989ca59a07afb535b0a3f","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"Northern Colorado Land Use Changes","cells":[{"id":0,"value":"$33","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":584,"value":"viewof cameronPeak = Inputs.button(\"Zoom to 2020 Cameron Peak fire\", {reduce:()=>{map.setView( new L.LatLng(40.607,-105.463),11)}})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":607,"value":"viewof eastTsome = Inputs.button(\"Zoom to 2020 East Troublesome Fire\", {reduce:()=>{map.setView( new L.LatLng( 40.249136,-105.95558),11)}})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":612,"value":"viewof windsor = Inputs.button(\"Zoom to new construction near Windsor\", {reduce:()=>{map.setView( new L.LatLng( 40.467584300504214,-104.96183395385744),13)}})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":618,"value":"viewof windsor1 = Inputs.button(\"Zoom to Upper Yampa Reservoirs\", {reduce:()=>{map.setView( new L.LatLng( 40.0483142058281,-107.06378459930),15)}})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":645,"value":"

Map of Land Use Classification

","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":53,"value":"viewof year = Inputs.radio([\"2017\", \"2021\"], {label: \"Select Land Use Data to Display\", value: \"2021\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":21,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":610,"value":"

Sankey chart of Land Use Change

2017

2021

","pinned":false,"hidden":false,"mode":"html","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":380,"value":"chart = SankeyChart({links: Data.sankey},{\n nodeGroup: d => d.id.split(/\\W/)[0],\n width,\n format: (f => d => `${f(d)} Acres`)(d3.format(\",.1~f\")),\n\n // Had to fork the built in Sankey Library so I could make the colors match the Land-Use classes\n adjustNodeColor: (c,t)=>\n {\n\n switch (t.split(/[\\W\\n]/)[0])\n {\n case \"Water\":\n return \"#419bdf\";\n case \"Trees\":\n return \"#397d49\";\n case \"Rangeland\":\n return \"#e3e2c3\";\n case \"Crops\":\n return \"#e49635\";\n case \"Built\":\n return \"#c4281b\";\n case \"Bare\":\n return \"#a59b8f\";\n case \"Snow\":\n return \"#aaeafd\";\n case \"Flooded\":\n return \"#7a87c6\";\n }},\n height: 300})\n\n","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":365,"value":"## Supporting Code","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":250,"value":"$35","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":36,"value":"html` `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":18,"value":"L = require('leaflet@1.9.2');\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":255,"value":"// github sets the wrong mimetype on js files so we using jsdelivr as a proxy to load this\n// after updating github i have to hit https://purge.jsdelivr.net/gh/grahamsz/leaflet-image-layers/leaflet-image.js to clear the cdn\nLeafletImage = require(\"https://cdn.jsdelivr.net/gh/grahamsz/leaflet-image-layers/leaflet-image.js?2\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":393,"value":"import {SankeyChart} from \"88895dac303f3ed9\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":361,"value":"import {PieChart} from \"@d3/pie-chart\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":226,"value":"{\n \nmap.eachLayer(l=>{console.log(l.name); \n if (l.name==year) \n {l.setOpacity(0.7) } else if(l.name==\"base\") { l.setOpacity(1); } else\n {l.setOpacity(0);}});\n yield \"Map Layer Visibility Set\";\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":400,"value":"import {d3Sankey} from \"@d3/sankey\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":429,"value":"Data = performAnalysis();","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":443,"value":"//d3 = require(\"d3@5\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W36","stdlib":"1","titling":"auto"}} 32:"$31:body"