DataPhilly Jan 2021: Satellite Imagery Analysis with Python

great thanks uh thanks a lot robert i think that's the first time ever um described as a cracker jack data scientist i really like that um so uh i'll get started um share my screen here um all right so it looks like we are good um so a couple a couple uh notes before we get started um there's i think there's been a little a lot of talk before about um the requirements for following along with uh this workshop um uh uh i want to first say that uh this does require some dependencies that are a little bit finicky to install so i really uh apologize for any difficulties people had um at this point uh the the main way that i was anticipating people would run this would be through a docker environment which i set up the resources to create um that does of course require that you either have docker installed on your machine or you install it in preparation for this workshop so um it is a little bit tricky to install at this point um if you haven't already installed it um you're certainly welcome to try right now but i'm guessing that it will probably take more time than you have uh then you'll need before the workshop gets started so i'm just a reminder that you're welcome to just follow along i'm going to work through everything and the repository that has the notebook and the data will remain up so you're welcome to um install the dependencies and um and run through it uh and run through it at your own pace later on um if you have docker installed but you haven't already built the image i recommend that you do that right now the instructions are right here this docker compose build command um you might want to do that and that'll take a few minutes to run so maybe you can uh run that while i'm giving a little bit of a spiel and showing some slides and then you'll be ready to go when um when we start with the notebook um also i just want to say i'm i'm really happy to be here um i think that uh i first committed to to doing this workshop in um december of 2019 uh which if anyone can remember uh what the world was like back then um the plan was was for this to be in person and um and uh you know we'd all be in a different state of mind but uh after a lot of time in limbo i'm really excited to finally be doing it even if if we can't all be face to face in the same room but luckily that does um allow for some of the folks that are not in the philadelphia area to be able to join which is really great um i am joining from philadelphia um but i'm really excited that people can be from all around uh all right so i um going to do a little bit of an overview of some just big picture uh concepts and uh and then we'll get into the interactive portion um so this uh [Music] this workshop is is uh on working with uh satellite imagery and python uh we're gonna work with not only satellite imagery uh actually the the um imagery that we're working with is drone imagery but it functions the same way as as robert mentioned it sort of all under that umbrella of earth observation but we're going to work both with imagery and what's called vector data which is geospatial data that uh is used to annotate uh the globe and we'll work we'll sort of learn how those two things work in concert um go over a couple things first i'll tell you a little bit uh briefly about myself um uh and then a little bit of background on this um earth observation industry for lack of a better term and and um some of the reasons why you should care about this uh being able to do the the work that um i'm going to get into and then uh we'll do the workshop uh get into some basics and um a little bit more of an applied section on um basically uh working through the data and showing you how to sort of prepare data for a machine learning um project so basically um you know i think the objective of this workshop is uh it's really the um it's really the work that we're gonna do is really quite simple um i'm you don't need to know anything about uh geospatial data or about um satellite imagery or overhead imagery or anything i the the objective is just to sort of give you a sense of um some of the different um tools and how they interact and then you know even though the things that we're learning are pretty trivial i'm hoping that you can sort of imagine how uh knowing how to do these things would enable you to do to sort of expand the possibility of all the things you could do working with geospatial data and imagery uh but first a little bit about me um hey my name is simon i i work at xavier which you now um all know about uh it's a it's a it's a really terrific place to work i'm i'm really fortunate to have found a home there i've been at xavier for about uh three and a half years i'm a data scientist um i do uh i work on a variety of different projects some of them um are sort of more traditional statistical models that incorporate uh geospatial data and it's about health or housing or all sorts of different things some of them are um and then increasingly i've started to work more and more with um uh geospatial imagery which we as a company have um like a really um significant investment in um and some of that has been um on machine learning projects which we have developed some software called raster vision that um is a framework for performing computer vision tasks on large format overhead imagery um and again i've worked on a number of different projects so um yeah that's about that's it about me again i'm i'm coming from i'm coming to you from philadelphia before i worked at xavier i uh got a master's degree in city planning which is where um i first learned or sort of got exposure to working with um geospatial data geographic information systems which is sort of the the broad term for working with geospatial data and that led me to azavia and i've been able to really learn a lot and grow on the job at um at xavier so just a little bit about [Music] earth observation data um and this sort of industry you know i'm i'm um uh i know my little portion of this world and i'm no means an expert in uh in everything in the field but uh the basic gist of it is that uh recently it's become a lot cheaper to put satellites into space and there's been just like a boom in the amount of imagery that is available for analysis um and then as that's happened um you know deep learning has has had a bit of a boom time and machine learning and um basically the idea is that like we have all of this imagery uh that is coming off satellites that um is more than humans can comb through all at once um and we're able to so like we're able to to develop algorithms that can automatically detect different phenomena within the imagery and people are just sort of starting to scratch the surface of all the different possibilities of uh things that we can do with it so that's something that um we at azavia are really interested in and we're trying to figure it out along with everybody else and and we partner with different organizations that are either have imagery or have um have uh complicated problems that they want us to try and solve and it's it's um a lot of what we do i i would um direct you to the blog at azivia.com we've um myself and some of my colleagues have written up a lot about um are different musings about um what's possible and what we've been doing and sort of what we we see as the future of machine learning and uh geospatial imagery and along with all of this uh this data uh there's a really robust ecosystem of um software both proprietary and open source that we uh some of which you will learn about today um that we contribute to and um are really like uh interested in advancing the state of the art of um so there's a lot there's been a lot written on this topic um if you're interested in that i highly recommend going down a rabbit hole about um all that this um industry entails um it's really fascinating but uh i'll cap it at that for today and we'll go a little bit smaller board because um if you want to work with uh you know work in if you want to sort of jump into this field and all of the work that's going on with geospatial data and imagery it's really important to understand the basics of how to work with it so that's what we're going to focus on today [Music] but uh before we get into that uh there's like there are a number of different um free image resources there are uh a lot of uh certainly a lot of uh satellite imagery providers that um charge a lot of money for their data and and um that's a whole different world but there are places where if you want to if you're interested in what we do here you can um go download data for free and uh some of it comes with accompanying labels uh that you can analyze or start to work to start to do computer vision on um there's uh there's there's aerial imagery of the city of philadelphia on open data philly which i'm sure anyone who works with data um you know in and around philadelphia is is probably familiar with um there's also these a couple different um really large uh data sets that include labels uh of different um different polygons like uh buildings and ships and and cars etc um xview and space stand those are run by different organizations and have generally had computer vision competitions associated with them that have really large cash prizes those you know your mileage may vary on trying to compete with people who are doing the state-of-the-art in machine learning but they're really important um large data sets that you can work with landsat is is lower resolution satellite imagery but that's um from nasa and the usgs i believe um and that is um another free data source for um really large uh satellite imagery and then hand in hand with the imagery um uh we'll we'll learn about like sort of what i'm talking about when i say labels geospatial labels um open data philly uh has lots of geospatial data about um the you know the philadelphia city of philadelphia and the the region um and uh you know even if you're not working with imagery that's just a good place to to um to start to find some geospatial data that you can work with um even if it's sort of just mapping or or the sort of work that we're talking about today and openstreetmap which is uh actually where the labels that we're going to work with today come from and for those of you who aren't familiar openstreetmap is um is uh is it oh it's an open source uh free map of the world it's like the wikipedia of maps essentially so it has um data on uh you know buildings roads etc it's essentially like a base map of the entire world um that people from all over the globe contribute to lots of other places so um if you're interested in this i recommend checking out those links this these slides are up in the uh repo so feel free to check them out all right so um that's that um i'm gonna jump right into the workshop at this point so for those of you who um have got the the um docker installed and the repository uh the repository downloaded and you've navigated to that workshop to that uh directory um it's i tried to keep this uh as simple as uh possible once you had docker installed so you should be able to either um just run this runs just run this run script um or use the command that's in the um in the um in the readme in the repository and uh we should be able to um open up this notebook uh which is uh that we're gonna hold uh this notebook has all of the analysis for the evening or all of the um just all the work both some annotation and the actual code that we're gonna be running so this is your one stop shop for following along with the um with the uh with the work i guess um and and um i i might uh take a moment to pause if anyone has questions um i don't know if ron are you planning to monitor questions and just sort of shout them out when they have them or should i stop periodically to to check and see if people have any uh we will keep moderating simon and i think one question here that's coming up most importantly is getting to this jupiter notebook right so one thing you might want to just mention about is at the end of the container as to which link they need to do to come up here because some people are getting no service token or you know they're asking for the tokens there so yeah if you just point to this yeah um i yeah i would try this um this lower one uh the lower one is that is 127.00 are people are people having luck opening that one uh i think people are missing that part so just highlight the part that you need to click on that 127 that's the one um i actually yeah i'm not entirely sure uh the difference between these two but um hopefully that works for people i'll give i'll give folks a moment to look at this and while people are doing that is google earth engine also a free source for satellite imagery yeah google earth engine is a um is a um it's software for processing uh large imagery at scale um so it's not um uh in and of itself it's its main function is not to be a resource for procuring satellite imagery although i do think that there are some of the data sets which i may have mentioned which i mentioned in the slide show may be accessible through google earth engine um so that's that's a project from google um in which you can uh do some sophisticated parallel processing to process really large data sets uh on their servers and that is that is free i think there's like a sign up process um i don't have too much familiarity with google earth engine i have used it but uh and i'm not entirely sure where it is as a project right now but that's definitely something that's worth checking out i can type up here about that um google earth engine uh was originally started by google.org the philanthropic arm of google and was primarily aimed at supporting conservation and environmental type projects climate change things like that um it is not a source for raw imagery it's they've ingested and processed it and made it optimize it for doing the kinds of um processing you can do in their platform uh it's a pretty functional powerful platform it has both a javascript and a python api but it's also non-trivial there's a lot to learn there if you're looking for raw data data sources for satellite imagery one of our colleagues joe morrison wrote a blog about this last year that you could find on our on our company website but uh amazon web services has a pretty great open data repository that has a lot of earth tree landsat and um uh as well as the european unions sentinel satellites uh it's a huge trove of data that if you're comfortable working in the amazon environment is pretty powerful thank you robert great um i i uh i'll sort of refer to you when you think um if you're if you're seeing the chat and seeing that there's um you know uh people are starting to get ready to run yeah yeah i i think you can start yeah yeah okay great are those people able to catch up we don't want to delay you sure sure thing yeah um great so uh again the um objective of this workshop is to like really uh learn the basics of um geospatial data how it's encoded and um how we work with it in python um or at least some of the tools that are available this is um i don't claim this is an exhaustive list of all the tools to to work with geospatial data and python but the ones that i've found to be really useful um so uh you know i i scoped this down to uh to make it as sort of simple and easy to use as possible we're gonna work with um two really uh just like two very small data sets that were um that should be quick and easy to to work with and we're small enough to just upload to github um a portion of a geotiff which is a a large uh image essentially with geospatial information um over um the city of uganda and um some a geojson file which is a geospatial data file um that includes uh labels of build buildings uh over that same area and those came directly from openstreetmap which as i said is the um this open large publicly editable open data set of all sorts of um information all around the globe um so uh you know i think that the point is to show you um just get you familiar with a few different tools um and uh get a sense of like all of their capabilities and i'm hoping that at the end of it you will be able to say um that you have a sense of the uh sort of possibilities that are available using these different data structures um and you could extrapolate those to do something um to do any number of things that are uh much more sophisticated than what we were able to get to today so um that's where i'm trying to get to and if that's uh if you have questions that you think will help um enable your understanding of them uh please feel free to ask um at any point and um i think just in the chat and ron will shout those out to me um so uh some of those tools uh specifically that we're gonna talk about um geopandas uh for so for uh working with vector data um which i will uh get into the difference between vector and raster data um at a basic level um geopandas which is a um which is built on top of pandas anyone who um works with uh data in python is probably familiar with pandas it's a um it's a library for working with structured data in data frames and it it provides us with these uh geodata frames which are um allow you to do all of the um all of the work that you do in pandas but also um enables you to do geospatial things um and then geopandas is built on top of a library called shapely which is essentially [Music] just a library that encodes geometries and basically allows you to uh perform operations on them whether it's finding the area of a polygon or interacting that polygon with another polygon or doing all sorts of different things that we're gonna get into and then um on the raster side which includes imagery um you know the two main tools that we'll be using are gdal which stands for uh geospatial data abstraction library which is a low-level library and command line tool for doing a lot of the similar uh similar operations with imagery and raster data sets it's totally indispensable um if you work with imagery you'll use it every single day um and then rasteria which is a a python library that's built on top of gdaw and essentially um essentially just allows you to ingest [Music] ingest raster data sources such as images into python and then convert them into numpy arrays that's the other uh really the third tool that i could have uh added to this list again anyone who's who's um familiar with working with data in python is uh is probably familiar with numpy which just uh allows you to do when you're algebra and work with like n-dimensional arrays of uh numbers just uh layers and layers of numbers which is really just what images are um uh so we're gonna do just sort of like input and output and then um windowed reading and writing and basically that's just like uh how do you read a portion of an image which is which is important because these data sets can get really large um and then in the second portion we're going to like apply some of the things that we learned in the basics to um basically uh turn these this geojson labels and this geotiff into a format that's like easily digestible by um a computer vision library so oftentimes if you're using an off-the-shelf um tool like fastai you know it accepts uh you know png data png images and like a csv of of like classifications for them or something like that so basically like how do you get from rather than uh the the way that you know geospatial data comes out off of um you know one of these sources so how do we get from from the input to to that output so um there's there are two different um tutor files uh scene.tiff that's this aerial image uh and buildings.geojson these are both uh in the repository so if you uh download the whole thing you will have those um and we're just gonna read and write directly from um from this repo um so before we go ahead and just um import some libraries again you have all these in the docker environment and then uh that will use throughout the throughout the workshop so um i guess the first thing i want to talk about is uh vector data so vector data is um uh in the gis world is usually referred to as points lines and polygons so it's essentially just um uh a uh it is is like uh data points that have um either uh you know area like a polygon that is like a closed loop of vertices um length in which it's like a line that has uh like sort of length but no area and that would be a line which which is often like how a road is stored as and or a point which is just like a pair of coordinates so um this screenshot shows you um what uh the uh a polygon looks like um in this format called geojson which is what we're going to be working with it's a pretty popular um geospatial data format for vector data and this is uh this is from uh just this like online tool called geojson.io which allows you to make a polygon um and it shows you what that polygon looks like um as a geojson this is just an example this is a fair amount where i live and um uh geojson uh is really simple it's human readable um it is um it is uh it is a json compliant format so it's essentially just made up of this like type classification and then this list of features and each of these features is uh a polygon um so or in this case it's polygon um so it's essentially just like a two-dimensional list of vertices that make up the extent of this polygon and those vertices are represented in this case in longitude latitude which i'm sure everyone's familiar with and that's essentially a um uh just just show it tells you exactly where on the earth's surface these uh this polygon is um so uh that's what uh that's what it looks like um if we want to look at um this this is uh just gives us a uh like a look into the buildings.json file that we're going to be working with um this is just uh for those who aren't familiar um you can do you just run like a bash command uh in a jupyter notebook using this this exclamation point so it's basically just like showing us the the first 30 lines of this building statue jason um so uh one thing you'll notice is that um we have these properties here so this is basically like tabular data that is associated with each of these features um i kept two different categories uh the building condition and the material um because we'll work with that later the other thing you'll notice is that um the numbers don't uh are clearly in a on a different scale than the um then in the geojson.io example and basically like uh what that means is that this is projected data and it's projected in a specific um coordinate system um so basically that's just uh this gets into how we represent a three-dimensional surface uh in a two-dimensional um like on a two-dimensional map um and the key to this is that it's in a unit of measurement that makes sense rather than latitude and longitude this is this happens to be projected in meters so um those uh coordinates actually have some meaning if you were to measure the distance between those two different points that we see that's actually uh a distance in meters um so that will um that makes it easier to actually do operations that are um that uh require you to use um measurements of distance that that makes sense so um we are going to uh we're just gonna start by just reading in this uh geojson file again this is this is really simple it's it's it will look familiar if you're used to um pandas but um it's just this read file command so basically what we get um this uh this uh g is a is a what's called a geodata frame um so this is a pandas data frame it is um it's a like tile class of of pandas data frames and uh you can do like any pandas data frame thing that you would be able to do on a regular one uh for example like filtering just the uh buildings that are wood um or um you know anything that you can imagine um there's also like special methods like the plot method in this case uh will actually plot um all of the buildings uh so you can see them visualized um i will also sometimes reference this um this you won't have access this this is a gis software called qgis this is uh just in this for this has a lot of different capabilities for analysis but for our purposes i'm just going to be using it to visualize things so this is an example of um all the buildings as these building labels as they exist in the world um and uh you just see we've plotted them here um so the thing that makes uh a geodata frame special is uh this geometry column um which is just a panda series made up of uh shapely polygons which i referenced so if you look at the first value in um the geopandas um geometry column uh it gives us this uh this is a shapely object um which if you uh like just return it in in a jupiter notebook it will it will plot the the boundaries and you see this is just made up of coordinate vertices um and this has their shapely uh polygons have um all sorts of different uh properties that um you know like this gives us this is this building is 664 square meters um but we just happen to know that that this is um that this data set is projected in meters um you could get the centroid which is like the center point of the um of the building and you can like apply um these methods over the whole data set so you could you know with one line of code we could get we could go from this geodata frame of all of the buildings in this area to a geodata frame um of just a point so just like two geographic coordinates or sorry two projected coordinates of the centers of all those buildings that might be useful for something that you're trying to do uh all right so we'll put a pin in that uh we're going to use assignment yeah sure before you move on to the next part there are a couple of questions related to what you're discussing um so regarding the coordinate system here that you were referring to for the projections and all that so where is the origin for that where's the coordinate origin uh like the you mean um like where is zero zero in this yes yeah yeah it's it's um it's just in a specific pla like i don't i don't know off the top my head like so this is um uh you can see the chord the like with a reference system for this so like gee yeah so this is like this is a specific uh coordinate system it's epsg three two six three six so like this is like just these are like utm zones so you could oops you know this off this information just like is stored somewhere um so yeah i guess like maybe the origin point i think okay it looks like it's like right here uh in in the it's right in the middle of republic so it's it's specific to this specific um this specific projection which like we have determined is most appropriate for the area that we're working in um so essentially like um you have to just basically because like there's no there's no um like obvious projection like like you have to choose the projection that makes the most sense because all of them are doing um they're imperfect because as i said you're trying to represent a 3d surface in two dimensions so you're basically trying to choose the one that minimizes error the most um and depending on which uh which um which which projection you're using it'll have a different origin point okay okay so that segues i think you almost answered that question the next question is is there a concept of a universal coordinate system or are there different coordinate systems to be aware of yeah um well i mean you know the there are the geographic coordinate systems which is like long la like uh we see up here um and this is useful for doing um sorry what we see up here and that's useful for like uh placing a um placing a uh a point on like at like in anywhere in the world like that applies to the whole world but it doesn't make sense in terms of uh measuring distance and such because it's not in an interpretable unit so uh sorry i'm taking a long answer to this basically there are some like universal coordinate systems that are good for doing specific things like displaying a um displaying a a web map like there's something called web mercator um but if you are um doing um like sort of measurement and analysis it probably makes more sense to go to use a coordinate system that's specific to the location that you're working in okay so does this mean that there is always like based on the region of the world that you're working in you could find different coordinate systems and that's what you would use for your projection so you can choose that in a way yeah absolutely like this utm zone there's it spans the whole globe so uh like there's a different zone for different parts of the globe um so yeah i mean this is a really big topic like you could you could um there are uh i'm not sure hundreds thousands there are lots of projections that you could use for various purposes and then people have um that have different requirements but um you know i i i don't want to go down too much of a rabbit hole because because it does but but um you know the good thing to know about this is like this is useful because it has it's projected in a um a unit that is makes sense okay thank you yeah sure all right great so so that's uh that's uh vector data which i think is like fairly easy to understand um raster data is a little bit um a little bit more a little harder to maybe wrap your head around but it's essentially just like um made up of uh it's like uh breaking the earth down into a grid um and then assigning a value to each um to each pixel in that or each square in that grid so you know the the size of the grid cell can vary based on your data set um but um essentially it's like um a continuous surface of values rather than uh discrete uh discrete points like you see with vector data so an image is an example of that of raster data because it is uh it's a continuous surface of values right um like uh um and that uh we could see that um by looking at this um using gdal this library i told you about by looking at the image that we're working with um i think i think earlier on i said the file we're using is scene.tiff i think i changed the name to image so just regret scene image.tif is the is the um image that we're working with so basically this um this gdal info command just uh returns a bunch of information uh summary information about the image that we're working with um and just uh so for reference this is uh this is the image um it uh shows just a portion of of this city um and uh it's pretty familiar looking it's an rgb image so basically what that means is that um it is made up of three different bands uh just red green and blue that um all have values between zero and 255 um and those uh those three bands represent are stacked on top of each other and for each pixel in this uh this grid of pixels that make up the image uh there are three different values that um represent a specific color so this is how all um you know this is just how images work uh it'll be similar if you just read in um a normal png image that you uh pulled off of google um but this um this uh gdl information gdl info command uh gives us a lot of different information that we can sort of work through so basically we're seeing that this image is six thousand by six thousand that is the um if you that is the size of the image in pixels we also see some information about the projection which we just talked about uh both of these um the image and the labels are projected in the same coordinate system um which i just i made sure was the case to keep things simple but um it's important to note that if you want to allow these data sets to talk to each other um they need to be in the same um they'd be in the same coordinate system um and um that is not particularly difficult to it's not particularly difficult to convert one data set from one coordinate system to another you just need to know what data set it's you know what is like the sort of encoding of the data set that it's coming from and what is the encoding of the data set that it's going to we're not going to worry about that right now because they are in the same coordinate system but um you know this gives you um this gives you a lot of information um uh pixel size so this is um this is the uh like actual size in the projected unit of the um of the of each individual pixel so this is 0.03 times um a meter it comes down to like a couple i think like four centimeter resolution um and that just tells you like literally just this one block of color would go in like this is a pixel this is four centimeters right so um that gives you a sense of uh that's the resolution of the image um and uh yeah so that's that's uh there's a lot of information you'll see that there's a fourth band and that's just um this is the alpha band and this just um you'll see uh when um we print out the data that um that is that just tells us like uh that determines which pixels are no data which is a classification that essentially like this is not part of the image um because these grids are um rectangular but the area that we're looking at um is not necessarily um it's not necessarily all covered by the image where we get too much into this because i've cropped the image to a specific um square there isn't actually any no data in there but that's just something to know about um so uh okay uh simon i have one one question here for you based on all the details that you were showing that was being pulled out from that particular image right uh in the top uh like without specific to uganda itself so does this mean that the image has to have been tagged inside does it need to have some metadata or tagged or it can be on any tiff image or yeah yeah the the the jews the information is like built into the file um i think is that is that what you mean like the um you know we could erase this uh geospatial image from the metadata and and it would just be it would it would uh it would cease to be a geotiff it would just be an image and that would be the case if we if we just converted this into like a regular png um you'd lose all of this metadata um and it would just be like a large png that doesn't know where it is in the world but that's what makes this geotiff geotiff right is that it does have this geospatial information attached to it so yeah it is um it is built into the file okay great thank you sure all right so um the uh uh we're gonna use uh this uh library called resterio to to read in the images and store them in memory as nonpirated so this is how the syntax works rest area.open this is a dataset reader object and that uh src.read returns um the uh the numpy array so if we look at this um this is going to take a moment that's all it is that once we read it in it's just a numpy array so um you know you can see it's uh four bands by six thousand by six thousand so you know at this point you can sort of um i might be coming to focus that uh if you uh know how to work with numpy and you understand the like all of the sophisticated math that you can do with it um you can do all of that on top of like uh with images um you know and and um that might not you might be thinking like well what is multiplying an image times an image that doesn't make any sense um but there are um like all these different algorithms that allow you to do um band math on the different um images to sort of like pull out uh to to classify water or classify vegetation or do any sort of um analysis that you um that you would find interesting so that's kind of the core of like um of like how you work with these images if within in python um so you know you could see like this is the the max value among this whole image is 255 again the values span from 0 to 255. um and then we can we can use this um this matplotlib function to uh plot the image in the jupiter notebook um so i guess this is flipped because i transposed it but um it's basically just plotting it on a cartesian plane um the six thousand by six thousand it's just um it just naturally this function knows how to plot a three band image um so you can look at it so um the uh that was easy right but you know you saw that it took um not that long but it took a moment to read in that image and um this is again this is just a really small uh portion of um an image that wasn't even that large so you know in the wild these things like these images can you know be tens of gigabytes um and it cannot be sometimes not practical to try to read the whole image into memory so that's where um this um concept of like windowed reading and writing comes into uh comes in and uh that'll be pretty key to what we're gonna work on today um so basically when we say like read a window of the image we're just gonna say all right instead of reading the whole thing just grab you know a one thousand by one thousand swatch from right here and let's just like look at that or you know let's let's do some analysis on that so um to do that you need uh you know we're gonna we're in that right now so basically um you do the same sort of there's the same sort of functionality where you um here we're going to open up the image and get some information about it um the width the height uh we happen to know these but i'm just going to store them as variables we know the width and the height are i think that's a typo there uh our um each 6000 that's a typo 2 uh rh 6000 um pixels i just fix that um so uh the um objective of this um and then this profile object is p uh just stores like a bunch of information about the information this is essentially like all of the geospatial information that makes it not just a numpy array um and uh right here we're gonna do a little bit of math to find um the oh right this is what we're trying to do yeah so we're trying to find the goal of this is to just find um is to take the uh the the 1000 by 1000 uh square that is at the very center of this um image so you can sort of see how this works right we take um the width of the image this incidentally because this is a square both of these values will be the same but we're taking the width of this image uh dividing by two that gives us the center point the center column in this uh raster grid and then um subtracting you know shifting that column over by half of the size of this window that we want so it's going to be you know sort of straddling the center of the image both horizontally and vertically so that's why we get these values which are just um are gonna be five uh sorry yeah 2500 right so um just just um smack dab in the middle um and the other thing that we need to do is um is uh update the geospatial information right because this um the uh right now this this like profile understands the image to be where it is in the world which is spanning over all six thousand of these pixels but like we want to write off a geotiff that has that has updated geospatial information so we wanted to know that it's just in the middle of this image um and we use this uh we do that by updating this transform object um so but first we need to um define the window and the window is essentially like um is this resterio object that uh basically specifies like where in the image we want to read from so um first off we're gonna um so basically we're gonna open um we're gonna open the image uh at this at the specific window that we've specified so essentially um this is just saying oh okay so um before we do that we uh are we we're to open the image and then def like define the um the uh carrot like define the transformation that will essentially give us the geospatial information associated with just the middle uh box that we're trying to take we're going to find that we're going to store that as x so essentially this is this this python object that just knows how to convert the geospatial information which is again is like where is this thing on earth uh from this thing is uh at the bounding box of this six thousand by six thousand image to this thing occupies just the bounding box of the 1000 by 1000 image at the middle so um we're going to put it all together we're going to read in um the this small portion of the image which again that this variable f uh will be a 1 000 by 1 000 band image um uh and then we update this pro this profile object um which is basically saying like okay you used to think that you were a six thousand by six thousand image that occupied uh this whole this large space in um in kapala you know now like understand that you are a 1 000 by 1 000 image that occupies a specific another a different specific geospatial location and then um we can you know specify this profile which says like okay now um understand that you're in a different place and then write off the smaller image so if we write that off you'll see that this um sample.tiff we just saved that as a geotiff in um in uh just right in this directory that we're working with and if we go to we go back to we can see you know you can see how it has essentially it's just saved um just this little portion that we want to work with and now we can do stuff with it um so that is uh [Music] you know the basics of windowed reading and writing so we're gonna build on all those things um in the next section but i guess i will take a pause i know i went through a bunch of stuff that if there and see if there are any questions um there's one question here i see what does each number and p transform mean yeah that's a good question i don't uh i've like never really uh committed that to i've never really learned that i assume that it's um yeah i actually i actually don't really know i've never really like needed to figure out what what each number um corresponds to but i assume that it's just like um some sort of offset math again this is just like a this is just a specific like python object that just like um just like it stores the information that you need to know but i guess this is the pixel size and then yeah it's probably it's probably like um like offset from origin in some way but yeah i've i've never needed to know that those are all the questions for now there were a couple others but robert took some time to answer in the chat oh cool great glad to hear it um all right great yeah just keep coming um if um if you if you have more um all right so like i said we're gonna use some of the some of this stuff um to uh just complete a couple pretty basic uh again pretty basic tasks um the first one we're gonna um we're going to work on is uh just like prepping uh some data for building chip classification so if you think of this um you know one of the three big uh beautification um tasks that uh that we do is chip classification so that's essentially um where we break um break an image into a grid uh and um train uh model to identify whether or not um a specific object uh there's a specific object present in that grid in that chip or grid cell so in this case we are going to um prepare for a computer vision task in which you want to um identify which portions of our image um include buildings you know if you remember um buildings are not you know there's some of the image has buildings in it and some of it doesn't so you essentially be doing like a building detector that that um gives you the general location of buildings if not um the specific outline or it doesn't return a bounding box the way an object detection problem would um so uh let's think about what we need to do for that um you know i know like a pretty simple way that um a computer vision library might want that is is basically like um two different um you know it probably expects uh png images of a certain size so just basically like the type of image you would download from google um or you'd get in a like basic uh computer vision uh you know training data set for for for like a simple problem um so we need to uh go through and break up this large image into a bunch of really small pngs and uh use the building label data set to identify which particular locations have a building uh covering a portion of the image and which don't and then filter those images into two different folders um so the way we're gonna do that uh first we're gonna read in the buildings file um this is uh the same geodata frame that we were working with before um and then uh oh yeah the first thing we need to do is um you'll note that the buildings uh cover a much larger extent than uh the image so the first thing we're gonna do is basically clip this building's [Music] data set to the extent of the image that we're interested in and that's a pretty basic gis operation you could think of it essentially it's just like a cookie cutter that you're putting down on top of the building's file and extracting only the ones we're interested in but um to do that we need to know like what how do we uh define the area that we're interested in um and to do that we need to find the uh bounding box of the image so the bounding box this is another bit of data that's stored in the image and we can open it from this dataset reader so you can see this bounding box um this i assume this this probably makes sense to you like the left most coordinate is you know 45 9 286 bottom 27 87 1 etc and um you could create a polygon from that using um by just combining the different vertices it's just like a pretty basic you can just think of like a square and a cartesian plot um but we'll extract this um this uh object into the different coordinates for left bottom right and top and now we're going to use those coordinates to create a shapely polygon to do some of the geospatial operations that we were talking about before so um you can see [Music] this is how you construct one of the ways you construct this shapely polygon objects which is uh just a set of coordinates left bottom left top so you're essentially just right top right bottom so you're essentially just tracing around the outside of the box and what that gives you is this um shapely object um and now um we are gonna take this shapely object and turn this into um a geodata frame so um what we get here is you know it's the same it's the same class as the um vbox is the um buildings geodata frame uh but it only has uh it just has like one row and it has no attributes other than just the geometry it's just it's just um a shapely geometry uh forced into a geode data frame so it's easier for it to talk to the buildings um and then you can see what that looks like plotted the same way um and that's again that's the old buildings and then we can use this um this gets into like these are some of the like special geopandas things that um differentiate a geodata frame from a normal pandas data frame and this is just like um this overlay function that does this uh cookie cutter operation that i was telling you about so now um this basically finds the intersection of um all the buildings and the bounty box so it's basically like um uh you just it's just like the place where both of those um both of those poly like data sets exist so now you can see the buildings oh if you plot it see it's just now it's like just the buildings that are um present in this um present in this bounty box um so right here um we're gonna create like a file structure that has like buildings true and false so uh i just added this in to just remove an existing um file structure if um you know you'd already run through this notebook so basically this like will remove the buildings directory if it already exists and then make um this uh structure that is like has a buildings folder and then a true folder and a false folder um we're gonna do that and then um we're gonna take uh um all of these buildings and um like we're not interested for the purposes of this task we're not interested in like each of the individual buildings we don't need to save each um png off with like the name of the building we just want to know whether any building intersects with that uh specific chip um so we're gonna convert um all of these buildings into one uh all of these individual polygons that you see right here into one multipolygon which is another shapely object that um just contains like a bunch of polygons so um that just you just take the series and cast that as a multi polygon if you look at that you'll see it's just like i guess it plots it in a different color but it's essentially the same it has all the same for the most part the same properties as the polygon and then um we are eventually going to iterate over the different um iterate over the uh the image in um in specific windows but um right here i'm just going to specify the um png size so let's just say for whatever reason whether it's like uh what you've determined to be uh the right size for your task or or it's like the um the computation the the image size that the library they're using expects we're gonna specify that as 600 and for our purposes that's nice because it fits uh neatly into six thousand so um these uh pngs that we're creating are gonna be um 600 by 600 um and uh so if you do uh six thousand by six thousand they're going to be a hundred of them right um because uh it's gonna be like a ten by ten there's six hundred uh a 600 uh yeah you'll fit 10 uh chips across the image across the 6 000 pixel image um so basically like right now we're gonna do um is uh go through and um you know which is sort of basically go over the image uh 600 pixels at a time read in each um read in each portion and then uh you know i'll actually i'll run it here because it'll take a second and we can sort of talk through it as we go so um basically uh so we're opening up the image we are uh finding the width and height which again are both 6 000 and then we're just gonna like iterate over that grid uh both over the width and the height um and uh and read like read and write one um a chip or you know sort of section of the image at a time so uh this is you'll remember this from the previous section uh we're gonna construct a window that you know at the beginning we'll start at zero zero and then be uh you know it'd be like 600 pixels tall and 600 pixels wide um but then the next time around will be at 0 600 and then span for the next 600 pixels um and then we are going to uh find the corresponding spatial coordinates because we need to know um like where in um like where in the world this chip is so that we can compare it to the building's polygon and determine whether or not they're like it is covered by buildings so to do that we need this transform that we talked about then we need to read in this portion of the image as a numpy array um so this number eventually we're going to save off as a png but first we need to perf like do that analysis to determine whether or not it overlaps with buildings so to do that we're going to do the transformation that we're familiar with from before and um get the uh bounding box of um the uh of this specific like chip um this this um you'll see here that we have to like open a temporary geotiff even though we're not gonna write to it just to find the bounds this is uh i did some digging on this it's it's like a weird limitation of rust stereo that you can't just get um the bounding coordinates directly from the transform but um this is just like a workaround for now and then we will use um the same uh this polygon function um this is a different way of constructing a polygon from bounding coordinates um i think it just needs the coordinates being a different um sequence than um when you get them from geopandas but essentially all you need to know is that like we've you know gathered the information about where this chip is we've created another polygon and now we have two different shapely polygons that can talk to each other so we'll use this intersex method to check whether the bounty box of the chip intersects with the multi-polygon which includes all of the buildings and specify that as this has building variable and then use that to uh specify like uh the the label which is gonna determine whether or not the the building has or whether or not this image should be put into the like buildings true uh folder or um the building's false folder so um uh so we use that to construct the file path and then we're gonna use um this um this library pill to save off the image now we're going to save this off without any spatial information because that's what the deep learning library wants it doesn't care about uh it just wants pngs um all that we need to know is whether or not it intersects buildings so we're just going to save off each of these images to their appropriate locations and then i just gave them a file name that corresponds to like where they are in the grid cell um in the uh we're like the width and the height of the starting point so then what we have is uh simon i think we have a question here sure yeah yeah so it's about the bounding box is it uh what do you do in the bounding boxes given in the utm coordinates when the bounty box yeah i mean that's in this case the bounty box being given in the utm coordinates is exactly what we want because the it all that matters is that it has to match with the um building labels which are also in the utm coordinates right so um at that point it's like it's abstracted from you know all that really matters is that they're like speaking the same language you know because if it says like i'm in this place in the utm projection then it just has to match with the buildings right so that's what i mean when i say like where it is in the globe um it's just like in a language that the other data set understands is what i'm saying yeah okay thank you sure any other questions yeah i think this may be the same thing about saying so this would be i think from one of the people again your image.tf and the building geojson should have the same projection system yeah yeah that's that's that's exactly right they always need to have the same projection system or you could you could do something where you reproject them on the fly if that wasn't possible but you know that's just getting into like sort of how you structure your workflow um but now you can see what we have if we just look at these folders you know we essentially have broken down all of these images do not have buildings and the images that are in the true folder do so we could you know take take a look at that this is the first or the second option from the building's false uh so you say this chip clearly doesn't have building in it and this is uh should be class this should have a building in it which it does so um now you could feed these into um uh machine learning like you could trade a model on this that uh just these are like appropriate labels now you know which ones have buildings which doesn't so it could it could then identify um which uh which parts of the image um do and don't have um like do and don't have buildings in it and then cobble them back together and that's the sort of like um converting the um converting the um that machine learning process to geospatial imagery is like is what we're sort of doing manually here but also what uh our library roster vision which i've mentioned a couple times that's really what it's it uh it does uh that's the real value out of it so um so that's uh that's one example um the second thing the second case that we're going to consider um is uh if you wanted to do uh like a semantic segmentation uh test which is another computer vision task and that is different in that um with that um in that case you want you we wouldn't just be interested in specifying like whether or not this chip has a building in it we want to know at a pixel by pixel level whether or not this is building or not so um that is uh that's a different task it takes a different out like a different input and so specifically what we need to do is to create a um a raster surface that contains the ground truth labels of where the buildings are and where they are um so this gets into like um how we think of like raster surfaces are not necessarily just images although images are one type of them they really can contain any arbitrary characteristic a common one you see often is like land use classification so you'll have um a map of a certain area and each area has a value tied to it like each pixel has a value that corresponds to a specific land use like um you know woodlands or concrete impervious etc and um so it's that's a that's like a categorical value that's just encoded with an integer um and that's essentially what we want to do here at the end of it we want to have um an image that uh has uh uh actually in the case that we're gonna do we're gonna do zero values for like zero one and two which is um uh building condition is good build the condition uh bad i think or uh sort of unsatisfactory and then no building and i'll i'll show you how we do that um but uh that's the end goal and um and just uh because again we're gonna do that uh for just like a portion of the image again um but first uh before we get into that uh i have a function here and this just repeats a snippet of code from for up here which is just that workaround that i was talking about this function just uh takes um it takes the data set reader which is essentially like the interface to the to the geotiff and the window that we're interested in and returns the bounty box um as a shapely polygon so we're familiar with this code already um we're going to use that later so um i guess yeah so the the in this case um the uh we're gonna want to uh rasterize that's what this process is called this process of um like burning values into um a um a raster file at the specific locations of um like where that like correspond to the values in a vector file uh we're gonna rasterize just like a one thousand by one thousand section uh right at the top right corner of our top left corner of the image um so again um you know i think at this point you're really enough to see like how this window corresponds to the um like the that top left corner of the image um it's starting at like origin like the the offset for the rows and columns is both zero and then it's um it's going to extend a thousand pixels up in a thousand per thousand pixels to the right and a thousand pixels down so um first we're going to uh we're going to read in that portion of the image as numpy array um and then we're going to use this function to find the bounding box of this section of the image as a shapely polygon and we're going to also grab this transformation which again gives us um uh allows us to update the geospatial information from the whole image to just the small portion so simon one more question so this building data that you have uh did it was it pre-processed and labeled beforehand yes well the billing data came from openstreetmap um so it's extracted as part of a um this competition and um uh yeah i did i did some basic i cleaned up some names and like variable names and stuff but um yeah it just uh it's it's it's straight out of open stream now okay um thank you yeah yeah no problem so again like uh to create this raster we need to know um like which uh we need to extract just the um buildings that we want to burn into the raster file and we're going to do that familiar operation in which we clip one [Music] clip like a geodata frame to a specific polygon um and actually i i showed you like kind of a more cumbersome way to do it above just because i thought it'd be useful to understand how you could basically turn a shapely polygon into a geodata frame but this is actually an easier way to do it um this geopandas clip function which essentially this rather than needing uh two different uh geodata frames it can take just a geodata frame and a shapely polygon which again just keep in mind that like this is all um there's not like um these are related right because uh geopanda's geodata frame is just uh um it's just the thing that makes it a geodata frame is just a bunch of shapely polygons so there's not a whole lot of abstraction there um so simon one more question sorry to interrupt you right there because we are trying to follow you through like running it on our own uh set too on our own uh notebook here sure these buildings that you have here right it goes back to a specific cell and i just want to check your gpd.overlay it is a little higher above sure yeah and it's not i don't know if anybody else is having problems but when i'm running it it is not working the gpd dot yeah because yeah it gives the problem about uh spatial indexes required either archery or the pi geos are you running in the docker environment yep in the docker environment and uh i don't want to take time away from you but that was a very quick fix so we could follow along i'm kind of uh uh it's possible that um uh rob is this is this um have you updated like did you build the docker environment like several days ago uh yeah i think i think i've actually i think i might have actually updated it since then okay then fine then let's not take time away from you yeah yeah so if that's happening i would um if you built it a couple days ago i i tried not to do this but i did add a couple things to the docker environment so if you just rerun the build script that should fix it um uh yeah sorry about that i know it is okay yeah um so uh where were we uh yeah so this clip is just like another way to do this we're going to do this like cookie cutter operation to extract just the buildings that we're interested in that uh overlap with the bounty box of the window that we want to work with so that does the same thing so you can see uh this is it this is just the small the small portion of the image these are the buildings um so if we look at these if we look at these um buildings we could see uh we see that there's actually a variety in the uh the condition that we're working with so let's say this uh a basic example of um this task might be to um predict whether or not there is a building in this um in each individual pixel but let's say it's like we want to do something a little more specific sophisticated where we want to classify the condition of the building based on the image now this may or may not uh be possible from the image because the condition may or may not be uh there might be enough nuance in that you may it may just be something that you can't see from above but let's just say that's something that we're trying to where we want to experiment with um so what we need to do is like burn that specific value into um into the raster data set um and um you can't use like a string value so we want to specify another um column that just like uses integer values to um uses like integer values to correspond to each individual condition so then we can do like a a familiar pandas operation here that just uh uh that just allows us to create a new variable and we're going to use this variable as the value that um specifies like good versus poor that we will like uh feed into our machine learning library in the encoding of values so um the way that we're like what we're gonna do to output this um this final raster is to use this rest stereo function called rasterize which does um which does this process that we um that i've been talking about so um you know i just uh like we'll show you how we just build up some of these parameters um so what this function accepts um the first thing is the shapes which is like um a [Music] um it's a like a a interval like a list or a um just something you get it right over uh that includes tuples of pairs of the geometry so like the shape of the building and the value that you want to burn into the raster um so uh i'll just do uh some i'll just sort of show you how we do that um these actually i think these are these functions are are um are not actually relevant i wrote these before i realized that we didn't need them but you just see where we build up the shapes here so we're essentially just combining um the uh this allows this is a another thing that's familiar from pandas where you just like combine two different series is um i guess this isn't even that's really specific but like we're just creating um a list of uh tuples that's um a combination of a building geometry and the condition value that corresponds to it so each one is this is the boundary of a building and the condition two um which is uh the um which is the the value that corresponds to uh a good condition building um and then we uh the second uh i think we pass in is this uh the shape of the array the shape of the array that we're gonna write off as um as a raster um so uh basically like you keep in mind that this like uh turning the uh the buildings into a grid it doesn't have any like inherent understanding of like the resolution you know you could you could make there's like there's no limit on how fine you want the rest you could make the resolution to be because these um these uh lines have no width but we want the resolution to match the image so um it needs to understand like okay how how big uh like what should the resolution of the array that we're we're rasterizing into b so we want that to match the uh the image the the section of the image that we read in so that's a thousand by a thousand um and then this fill value is the value that we want to fill all parts of the um of the window that don't have a building on top of it so we're going to specify that to be zero the transform is uh what is going to um it's gonna allow us to convert the geospatial information to match the specific window that we're looking at and um so i should actually run these uh and uh this is just uh the data type uh we wanted to match the image as well and then we can run um this function with all of the different parameters that's uh [Music] oh whoops so um so that returns uh a um another numpy array which now has the values burned into it so you can see i just found this center of the image you can see this gives a little summary but she's like this is showing different portions of the image that have no building have a building in poor condition uh building in good condition and that's like the sort of the edge of the building that is in good condition bordering on the area that has no building so we can look at that so again like the the difference between this and the image is that this is only one band it's not uh three different values that combined to return a color it's just showing the labels as a flat uh single band raster um but if we could still look at it and here i think it just chooses these colors arbitrarily the plotting library does but we have the different places um you know if you were to zoom in you see this is like a raster grid of pixels um and the different values are represented by um by uh different colors so if we look at the image to match it you could see how it lights up right these buildings with the red roofs are the buildings that are in good condition this one right here with the gray roof is in poor condition and then everything else is non-building area um so then um you know you could just write that off as a geotiff um or as a png depending i didn't go through that step because uh we've sort of already learned how to do that um we're getting short of time um i can briefly show i'll briefly just show this this final little uh step this just shows how you could do the same operation um just using a gdal command um so uh but before we need before we do that we would need to take uh the buildings that we're interested in rasterizing which is um just uh again we've specified those buildings that like overlap with the the small window and um and also importantly we've specified we've added the condition the numeric condition column um so we're just going to save that geodata frame and this is this is how you save off a geodata frame to a geojson or another spatial file like file type that um a software like this qgis that i'm using can ingest we're going to save that off as rasterized buildings.geojson so for the buildings the buildings that we're hoping to rasterize so you can see you know if we bring that in you'll see it's just it's just the buildings that we're interested in that just cover this specific location and then this just shows how you can do this with um a genoa command from the command line so this is going to do exactly the same thing it's essentially passing the same information it's saying but it does a lot of the um it's doing a lot of the work for you because it sort of it automatically reads in the spatial information and converts it um so this is basically saying like rasterize um this geojson file at this 1000 by 1000 resolution uh use this condition num variable know that it's in this specific utm zone and then write it out to this buildingraster.tiff file [Music] so do that it gives you this little progress bar and then um you can see right here that this is um you now have uh this raster file that's just represented as zero zeros ones and twos so now you can put this into and again it's this is pixels right this is not the um the vector label as we see it um before it's rasterized so this is now in a format that you could use for um a semantic segmentation test like i described before um yeah and then i guess uh oh i guess i ended up writing off the thing from the step four um but yeah you could see um oh sorry so this is just yeah so this is just reading in the output of that file and you can see that it matches matches up with this one um so that's uh that's it for the workshop uh i i think we've got a little more time time i think we we said we're gonna go till eight so i'm having to answer any more questions that folks have um but otherwise i really hope that this was um this was informative for people thank you very much simon this was really good and helpful so if anybody has any other questions you could unmute and you could ask simon directly if you want

Channel: DataPhilly

Views: 2,892

Keywords:

Id: KFO9zCFBwtY

Length: 98min 27sec (5907 seconds)

Published: Fri Jan 22 2021