Analyze floods using ONLY Python! (aka spatial data science)

everyone so today we are going to go through a quick tutorial um creating a really cool map in Python it's called a River elevation model and basically what a River elevation model does is it takes the height of the river or where the actual kind of water level is and then shows all the heights above that essentially what it does is it gives a very clear picture of all of the movements of the river and the very slight changes in elevation that take place from the riverbed and you can see how Rivers change and flow and move over time and depending on the river if it's sort of in a very flat area you can really see this and there's a really good example of this from classic cartography you can see this map here which is a lot of people know this and it's pretty famous for showing this someone actually went out and mapped all this to show how the Mississippi River changes and moves and and you know goes different directions over time so that's what we're going to do today I'll walk through how we kind of got there we're gonna be looking at a few different packages today working with raster data so we're going to pull a digital elevation model from the U.S Geological Survey pull that in and create our new map from that so let's go ahead and get started so all of this got started kind of from a post I put on LinkedIn I had you know explored this package a little bit and really anything related to hydrology Rivers um you know even raster data is not my strong suit I don't do a lot with this and haven't so this is really kind of new learning for me but I found this really amazing package called High River and this visualization like really struck me if you see this in the example so there's a lot more you can do with high river this is not all the things that it can do there's a lot here this is a python package of course so you can look at you know flood station samples you can look at land use and sort of the Watershed of an area um you know dams splitting Ketchum areas there's there's a ton you can do here your River elevation across section for you know hydrological use cases both rivers and you know coasts and I mean literally everything it's a very very cool package and really contained but I mean the first thing you land on is this awesome map um and and you know you can certainly do this with this package I it's a very you know similar type of use case and you can plot that here and you can see it looks really great um what I did find was actually this you know post over here which is actually sort of another um another map you can see there's some different things here there's a tutorial if you open that up you can see exactly how they did that here very similar you can see you know there's a whole story map that walks through how you can do some of these different things that one's gonna take a while to load um and then this is The Notebook I found which I have open in another tab here that I actually started to pull and play around with so a lot of the thank you goes to Daniel John if you check them out here uh you know I actually love this description on this GitHub landing page so you know uh Daniel put together this tutorial Daniel put together this tutorial I I think it's amazing he does a lot and posts some different things um about you know digital elevation models and geospatial in general so definitely take a look at what he's up to um but going back to the tutorial this is sort of the same type of uh process you know relative elevation model and I basically pulled off of his notebook to build the you know the tutorial we're going to go through today which is all in this notebook here so um I can I'll put this on GitHub so you can share it I'll clean up the descriptions and make it a little bit easier to use um and then we'll walk through that but before we do that I'll just kind of talk through a few different libraries that we're going to be using here um just to get everyone familiar so one is rest area we're not like directly using this Library a lot here to be honest um and we'll get to that in a minute but rest area's ability to read open and do some visualization of raster data so there's a lot you can kind of take a look at here you can read a raster file well you can extract the different bands visualizing it do some interpolation so lots of cool stuff you can do with that the other is is x-ray which is not a geospatial specific Library it's basically reading python you know lists or you know arrays if you take a look at the home page here this really makes working with like multi-dimensional data in Python really easy and of course rasters at the end of the day once red and python basically show that you might have a you know a list for a specific band that had an index point and then you might have that same index over time right if you want to detect time change so it's become pretty popular to do a lot with that and then there's this kind of Blended version of this which is uh Rio x-ray which is basically restorio plus X-ray and some basic functions to pull in some of the functionality from oursterio like reading data here and then you know other things of clipping or reading you know Cloud optimized geotiffs and lots of other stuff that you might use x-ray for so it's kind of this in the middle Library so a lot of the stuff you know is pretty standard numpy geopandas pandas osmnx osm and X if you haven't used that before we'll use it for one very specific thing uh shapely sci-fi spatial which basically has our um you know inverse dataway calculations through KD trees and then some plotting data Shader in that plot lib so that's what we're using today to do all of this so let's go ahead and jump in um the you know if you haven't loaded your libraries yet you'll have to install them you can actually install them right in the notebook here I recommend using this Dash queue after your pip installs it's called it's like the quiet flag so you don't get all of the readout basically make your notebook super long when you install it if you haven't done that you can go ahead and just hit run all your libraries will be imported and then you're ready to start going so that's to get started so the first step is actually to load a digital elevation model and there's the USB GS data Explorer for lidar and you know in the you know Spirit of kind of creating a similar map to um what we had you know the original Mississippi River map direction and take some data from there um so if you what all you have to do here and the area that we're going to sort of map is around this city cape giordo in Missouri right it's right on the Missouri uh Illinois border so it's kind of a right just before the Ohio River over here joins up and flows into the Mississippi River at Cairo uh Illinois right this tip here so there's some interesting change patterns over time and then you also have these Hills here with these national forests that kind of Where the River clearly cuts through and then opens up into this plane down here and then so let's just go up here so it's super easy to use this application we're going to use a Dem and you just Define an area of Interest so all I have to do is so all I have to do is just kind of hold and drag an area here and you'll get a bunch of the Dems there's lidar data but then you also get the digital elevation model data within that location and I used I end up using the 10 meter data you can download those here and if you just note they're separated by dates this is uh you know it looks like March 18th 2022 uh you know August uh or July 21st 2022 so on and so forth so you just pick the date you want and then you can pull those and they have 30 meters as well and they have very you know small very you know one meter level Dems as well so these I'll show you why I ended up working with the 30 meters because you know there's some of these issues where the scans aren't complete right so you just get all these zero values in the side we'll go to maybe one near a river here it's kind of annoying you can't actually click on the map to get what you want so you see all these all these zero values and then you just kind of get some of it in some of it out so that the one meter I kind of said okay forget about that similarly with the and it's similar with the three meter two it's the same type of thing so we can pick this area here it certainly works but then you just get these vast empty slots depending on the you know the time and you want but the 30 meters are pretty consistent and we'll just let's for example take a preview look at this one here and you can see yeah this is pretty much the same coverage now if you want to use a different file I also was testing out some data from uh not the USGS but another state if you get these extreme zero values and I'll actually show you what that looks like here in the first step um so let's take a look at this looks like in the first step I have this data in Minneapolis so we'll just copy this out for now and we'll go with the as one meter Dot and if you see this and we just try to basically show what this looks like so if I load that in I'm gonna this basically shrinks the Dem and you can kind of see here I'm just going to get these this massive set of like zero values so nothing shows up and and what I did was pull that into qgis I saw where all the values were and then I ended up clipping that but I think even then these Dems that we have from the USGS cover pretty large area right we're talking from you know uh I think this is like a full degree of latitude or longitude something like that so you got to be aware of how big these are and pick a area to focus on otherwise you know there's such extreme values in some of these Dems so I actually recommend clipping this anyway and what we'll do here is do that but with our original value so here we go we're going to load our geotiff that we are using this is the uh you know this is from 2023 in January 17th you can see here it's North 38 West 90 so that's your kind of you know where this fits on the grid of the Earth as it were um now what I did is actually clip that so I'm going to import geojson you might have to install that Library as well and then all I had to do to get my coordinates was go to this wonderful tool if you haven't seen this before it's just bonniebox at cloakingtontech.com I use this for just about everything you can search anything from osm and it basically just takes you over there and starts you with a basic bounding box so I kind of guessed at the area I wanted I don't remember exactly what I drew something like this and then you can pick your format like if you want this in just the you know the raw values in our case we want the gojson which you can copy here and then just load that into notebooks that represents our geometry and then what I did here was picking basically folding our geometry into an array right so we have our string geometry we load that as a proper geojson geometry here and then that represents our cropping geometry and then here we actually use that you know the combo resterio Library if we go back over here to Rio x-ray and we have these functions that come along with that so there's the Rio they call them real accessors basically here and there's one forklip and that allows you to clip the geometry by another geometry image that's all that that is we know that this is in the 4326 projection so run that and then that should just clip our Dem to the appropriate size and scale that we want all right so that once that's done this is basically going to run a process called a coarsening and trimming to kind of compress the Dem because it's a pretty large file at the end of the day so this does that um we're basically going to clip that there and then we just plot it to see what it altered looks like and here it is so once we have this plotted you can kind of see there's these Hills and it kind of becomes clear and you can just barely see where the river is which is why we're doing this whole process to get figure out what's in here and to really get that to pop out in the imagery so the next step is we need the coordinates of the river and this is where we use os7 and X which is how we've imported this library and there's this nice little function called geocode to gdf what I'm geocoding is the term Mississippi River um and that's all it is and then you can basically specify the result type I found that one work here you don't need to specify this and basically here I'm grabbing the coordinate reference system from R once again our raster data and applying that to the river so we pull that we find the Mississippi River and then yes we see that this is the wholen Mississippi River starting up in Minnesota going all the way down to Louisiana here so to find the bounding box of the roster you can call the roster access this Rio function then get the Bounce from that which returns our bounds here and then all I have to do is pass that into a variable and then pull the X-Men or Y Min x max y Max off of that so you can do that here and then we can actually clip the river to the bounds and then we have to grab that geometry there's only going to be one so it just pulls that whole thing in and call it the river geom and then we just take a look at it here this is just some basic shapely functions that we're you know pulling in here you can see what this all looks like and we can take a look at the whole thing here so if you want to just call the variable for River it's a nice data frame just as expected coming back and you can see all the different data about that here Mississippi River Waterway River and then what we care about of course is the geometry and then to get that actual value we want to find that at positions there we see our indexes Arrow we use the pandas function I look to find that and that's about it so we'll get rid of that for now and then what we're going to do is basically you know slice our bounding box and sort of crop it again to get to those values so we're just accepting the Y and the x-axis is here and then we're just going to basically create a plot this is just matplotlib nothing fancy here to plot our two variables on top of each other so we can see we have our nice river running through exactly where we'd expect it to be right so now we actually have to create the River elevation model so there's a few different steps here I pulled all this from the original notebook so nothing new that you're seeing here as I said I'm not an expert in this field so I'm just kind of going and following others instructions as we go along so what we're doing is we're going to sample the elevation along the river and what this is basically doing since we have our our data in arrays or lists in Python we can simply just pull those off right we're going to grab the X and the Y and pull those off of the geom if you remember this is going to pull those off of the river GM so we Define this function we basically turn them into an X and A Y and then on the raster variable we're going to do the exact same thing so what we do here this XR is how we imported x-ray so it's just a shorthand if you want to use it like that sometimes you see import pandas as PD that's the same thing we're doing there and then we're creating this function data array and we're basically pulling the values off of that so we're going to have two separate arrays one for your latitude and one for your longitude and that's all we're going to do there so that's basically all that's happening here if you take a look at the actual raster itself that we have you'll see that it it has those values in there so let's just quickly look at that as well it's called cropped and you'll see that this is x-ray data array and you can see there's three different parts here there's the band which is the elevation the Y and the X so there's a you know basically a list of those and you can take a look at each of these see all the data that's in the band here's some of the index you can look at X see what's in there it's all your coordinates and things like that so all that data is already contained within this right and you have all that data just in basically lists and that's all you're working with which makes it pretty fast it's a very base data type in Python so it makes things pretty quick overall so our final step here once we've gotten that is we have to pull in the interpolation function we can see here that that's from X-ray and we're going to pull that interpolator around the coordinates and then find the nearest Value to those other coordinates and what that does is basically that finds the elevation level at the river that is going to be stored in this sample variable here and then what we do here is basically prepare this in sort of a three-dimensional array which is numpy v-stack is from what I understand is basically a three-dimensional array we're pulling in the sampled coordinate values and then splitting that out into X and Y valuables and then just transposing that in effect so we'll take a break these down a little bit and see what happens at each step this peels off the X and the Y from that variable up here so the cropped coordinates and that's cut so this is pulling off the X and the Y from the original Dem and finding those and then what we're doing here is just interpolating those new values off of that and then we're sampling them so this last step is basically turning a flat in the race we're kind of just doing some data manipulation here and we can kind of take a look at what's happening in each step so let's see what C sample looks like this is pulling from the sample data and you can see this is just an array with your latitudes and longitudes along the river right so that's that first step in the next step here let's just drop this below so we don't have that same issue again we can look at C underscore X which is basically an array of your X values and then c y should be a array of the Y values and then finally we prep this data as you can see here for our interpolation so this reshapes the Dem we have our sample and then now we can start to interpolate those you know elevations from the Baseline of the river outward to the other elevations so when we're ready to do that we use that final function as the KD tree that's going to sample our sample the variables here and we can basically pull that up that's going to establish our KD tree which if we want to see what destruction python is you can always go like that and that's just a uh from scipy there's a set of spatial functions and that preps that for Kate so if you want to take a look at what that does it's a k dimensional points and you can basically it's kind of like an index to rapidly look up Neighbors and other things like that across you know different variables so this is where the interpolation happens you can basically pick your samples here so this is you're going to pull the tree you're going to query and interpolate off that interpolated values this is your sample which is the river once again that's the elevation of the river and then your interpolation which is everything else now you can kind of you know pick this or you know having some weights about the distances and adding some different you know basically math to to and you know calculate this now I'm like I said I'm not an expert on all this but um you know kind of understanding that there's some you know interpolation and some extra math going into this as it goes along to reshape the data to visualize so we can do here is we'll actually run that the original notebook starts at a measure of k equals five so let's just take a look at what that looks like we'll recalculate the interpolated values in the raster and what's really nice here is just you're just working on a raise and then you can render that basically as a raster and output that so you can see it's pretty sharp here from the interpolation about the relative elevation as you start to go up that'll get a a little bit you know um smoother I should say I mean you could go all the way up to 50 and see what happens here it's just going to kind of smooth everything out I found that the at least of the level of data that we're working with and maybe it's the size of the scale that worked a little bit better because we might have some Extremes in our elevation and if we know from the beginning we have some hills or mountainous areas that we want to render appropriately so that's what I did here you can take a look at all these individual steps if you want to see what the weights are you can kind of see that those are the weights along those array so with the weights we're just kind of creating the reciprocal here and then dividing that by the you know the ultimate total so it's sort of a even waiting scheme that we're creating here and we do that we take the weights times the values at the different levels and then that's what pops up here so you can see the actual interpolated values which is ultimately what we're going to map down here so we're creating a new data array from our interpolated values doing all that from the reshaping from the original X and Y variables here's our coordinates and then that's what we plot down here so if you're doing this and following along I would depending on the data you're looking at a play around with that K value and figure out what's going to work best for your data but this is sort of what I settled on and then this is the final step really simple we just basically subtract the original data the cropped data from the elevation raster which we've created here so that's all we do and once again it's super fast it's just lists of data it's just python list nothing else which is really cool and now's the fun part so this is the map that was originally created and what we're going to do this is without other elevation data in it you can change your color scheme if you want to here but we're just going to run this and here we get exactly what we wanted so this is our River elevation novel now it got flipped upside down at some point if you remember all the way at the top are the bottom of that Loop is actually South so you'd have to kind of rotate this to get the original values but you can kind of see where all these things you wouldn't notice in the original value were from water Flowage connecting here areas where you know the river used to run in these planes run off I mean these could have been floods it could have been yeah I don't know it's just very very cool and if you remember back to our original map this is kind of what was there so this is really cool to see just all of that highlighted and pulled out so the actual cities right over here somewhere I think if I'm remembering correctly and then the next last step is we can actually shade this along with the original Dem which gets this extra cool effect of the elevation plus the REM so that's all that this is doing we're using XR spatial and pulling in the hill shade to shade that and then basically stacking those with our original shading up here which is just creating it sort of like an alpha level or transparency so we'll run that and then that's where we get this really stunning map that I think is just awesome looking I mean you see these you know beautiful Hill Shades um you know we we have the river running you can even see down to like land plots here where there's you know farms and different things but it's pretty cool I I just I really love you know just looking at this and playing around with it you know this is something that I've never haven't really done a lot of but you know eager to explore and I think the underlying message here it's a pretty simple process to create these there's some other tutorials uh that Daniel calls out at the beginning doing this in huge is as well so you can actually do kind of the same thing here and this walks you through a full kind of history of Rems and what they are and you know uh floodplain visualization all this different stuff right so you know I definitely recommend reading this I'll put the link there as well so you can check it out but all right just so cool and so fast because at the core of it if you boil all this down with the data types we're using it's just lists in Python which make it super fast um you know it's a base data type and you can just do so much so quickly which is what I really enjoy so so here's the data flipped in the right direction after I actually did that correctly but I mean it's it's super the resolution is you know pretty good for this you can really see fields and things stand out as you start to explore this and then very clearly where the river's cutting and you know going into these Hills and these little streams that feed the Mississippi and then you know obviously where you know floods may take place and then these probably historical tracks of of where it's been before um I don't know just really cool I really enjoyed this you know kind of working through especially because I haven't done a lot with this before so that's kind of the tutorial um check it out and you know all you really need is to find a source of Dem data pull it in and play around with it make sure you're cropping it to the right selection but the rest of it you can kind of just run through you know you can dive into the specifics if you want to learn a little bit more about this but you know this is just a good starting tutorial that I had a lot of fun with so I hope you enjoyed this uh once again I'll share the notebooks and data sources and articles All In description so check it out there and yeah we'll see you next time

Channel: Matt Forrest

Views: 12,201

Keywords: gis, data science, spatial data science, spatial data analysis, geospatial, modern gis, matt forrest, esri, arcgis, gis technology, spatial analysis, geographic information systems, gis analyst, python, geospatial python, geospatial python tutorial, python gis, python course for beginners, sql tutorial, geographic information system, qgis, earth observation, hydrology, flood analysis, raster data model, digital elevation model

Id: EOzovrGOaJE

Length: 24min 6sec (1446 seconds)

Published: Tue May 23 2023