Processing Geodata using Python

Speaker 1 [00:04]

Thank you for the introduction. Thank you for attending. I will talk about geodata and Python. I will not show slides today. I'll go to the Jupyter lab I prepared here. So I think most of you already have to do something with geodata because geodata has spatial information and almost all data has some kind of spatial information. So what am I talking today? I do the foundation. If you want to do machine learning, deep learning, et cetera, with geodata, you have to know how to read this data, how to access this data. The big problem of geodata is it's almost every time very big. We have terabytes of orthophotos. We have terabytes of satellite imagery. And you have to know how to access that. And I will show some libraries today, some very important libraries. First, we have the GDAL-OGR, it's a C++ library with Python bindings. It's a quite old library created by Frank Wormadam maybe 20 years ago around. And the second library is GEOS, it's Geometry Engine Open Source. It's also written in C++. And as I said, there are some bindings, and documentation is C++ only. So it's great for Python, isn't it? No, it isn't. So some people created new versions, Pythonic versions out of it. And with GDAL, it's divided into raster data. So it means you have satellite imagery, you have images, and the second is vector data where you have lines, et cetera. and for the raster data we have rasterio it's created by mapbox it's an open source pythonic module which is really written in python so you have really python classes and so on and the second is fiona it's for vector data it's based on the ogr part of of the gdal library which is for vector data and also very pythonic it's also with classes etc And then the third is Shapely, which is based on Geos, which is also a Pythonic library. And now I move a little bit down here. You will see there is one very nice library called Geopandas, which is based on top of Fiona and Shapely. At the moment, we can't access raster data with Geopandas, but maybe that's something which will be possible in future. GeoPandas, something is missing, is built on top of Pandas. And then it's dependent on NumPy, et cetera, et cetera. It's just a short graphic here. So if something is missing, there are, of course, hundreds of libraries behind that. So we need Python 3.6 or higher for this notebook. This morning, Valerio Macho said something very important. We have to write the version of Python, which is required in a notebook, so I did that here You can try it with 3.7. It works. Actually. I'm running this on 3.7 Then actually I'm not expecting you to do that now. You can download this. I put a link here somewhere github.com martin-christen-pycon.de 2018 I will put that on Twitter too later if you miss it now remember it okay remember okay great so that's the installation you need this is just the top library so because there are dependencies like numpy it will automatically download this of course so let's start with shapely based on geos does one thing it can handle points line strings linear rings polygons multi line strings, multi-points, and multi-polygons. So if you don't know what this is, it's quite simple. Point, I think everyone knows what a point is. A line string is just a bunch of points connected, so you have lines. And there is a linear ring, which is just closed. I didn't put it here. And then, of course, you have polygons, and polygons are quite difficult because there There are some special cases with overlapping lines, et cetera, simple, non-simple polygon, convex, non-convex polygon, et cetera. And, of course, polygons can have holes. There are some countries which have holes, enclaves, for example, and there can be many, many such holes. And there are some rules. I'm not telling you the rules at the moment. For example, in a hole, you can't have one point on the edge of the polygon. In this case, you have to create two polygons. There are many, many rules. You can check this out. There's a standard. It's in the simple feature access standard from the OGC. It's the Open Geospatial Consortium, and there is also an ISO standard for that. So let's do something with Python. So you see three lines of code. first i import from shapely the polygon from the shapely geometry module polygon and then i just create a simple polygon we see here the coordinates and then you see that's why i said python 3.6 i make an f string i call the error and the length of this polygon and if we're lucky this will work of course it works so you see if you can access that directly so we can calculate some things but of course that's not all let's look at this polygon so Shapely directly has an output in JupyterLab or Jupyter Notebooks if you run this on Firefox it will not work because there is a bug everyone is pointing to the other because Firefox is quite strict and this is a SVG element which is not quite in the standard and the one says, okay, it's a bug in the matplotlib. The other says, okay, it's a bug in Shapely. And so we have to put two tickets on that. Then we create a second polygon. We can look at it again. We can create a union. We can create an intersection. We can really operate with that. And that's based on the Geos library. It's really a very cool library. Every big software packet in the industry, So if you buy a commercial GIS software, this is included, the GIS. So they use that. They sell it for big money. So symmetric difference, it's a cool thing too. And, of course, from such a result, you can calculate error lengths again. Or, and that's the more important part, you can create a thing which is called a well-known text, WKT, and it looks like this. It's a... Oh, I lost the connection. It doesn't matter. It's a multipolygon in this case. We look at it again here. You see this hole? Let me go down. And this one is text, and this WKT is a representation of a polygon in text. Text is very important because we can print it, we can store it, we can keep it forever. That's important with geodata. So, if you have an exotic binary format, we can't really store it for the next 2,000 years. So, this one we can write. There is also a WBT for binary, but I'm not going into that. So, we see, of course, that's string. And then we can load it again from the string. So, with Shapely WKT, we can load this and look at it again. We see, of course, it's the same. So this WKT is a very important standard. If you are using a GIS software like QGIS, then you can create a CSV with some information, and behind that you can also put these polygons, and then you can see it in your GIS, just by the way. There are some other operations you can do with shapes. You can check if a shape contains another shape. You can check if there is an intersection and so on I just made some examples here. It's it's a really cool operation those of you who Tried or actually did implement these operations. It's really it's really hard I remember many years ago when I was still studying we had to implement these things It's really it's it's hard to support all different kinds of polygons So let's go to the next library, Fiona, where we can actually read data. I downloaded here some examples. It's from the Natural Earth data set. It's open data. It's public domain, so that's why I put it also on my GitHub repository. It's no problem. And there I can just open this data set and look at the first entry. And we see it's called airports. So we see we have many metadata here. Here is the geographic information. That's the coordinate, the geographic coordinate of this airport abbreviation of the airport and so on. Even in different languages, you can look at that. It's a normal vector format. And then you can just access it usually like this array notation or, in this case, a dictionary. and access those things. So nothing really special, and don't forget to close it. It's important if you have Windows and you open too many files, it will crash horribly. So let me open this file and search for the Frankfurt airport. There are two ways to do that. The first one is I can just open the whole file. It's a couple megabytes, so it would be no problem. But let's think about the big data sets again. This could be 10 terabytes, and then you have a problem. So the way is just to iterate through the whole data set so you don't really need too much memory. So let me go to all airports and check when it's the Frankfurt airport. I will just print the property's name, geometry coordinates, and the property's Wikipedia. That's just a Wikipedia link to the Frankfurt Airport. And we see it's quite fast. I could even do a return here. So I could optimize it. It's not optimized at the moment. So we have the Frankfurt International Airport. And here is the coordinate. That's the longitude, latitude. And here is the link to the Wikipedia. I'm not opening it. I think you can guess it's true. Now I talked about geographic coordinates. every vector data set must have a coordinate reference system and I can just check it, what's this and we see it's EPSG4326 and EPSG means European Petroleum Survey Group they have a catalogue of the most important coordinate reference systems this one is a webpage created by someone okay, it seems I don't have internet at the moment which is not good for the later demos we will see about that so I connect again but you can check it out yourself you go to EPSGIO4326 it's called WGS84 it's a geographic coordinate reference system representing the Earth so internally it's using the pipe roach library where you can make transformations from one coordinate system to the other so let's look at different data set here i have the countries data set so every country of the world i have a geometry and some additional properties information so let me open it and look at the first entry again i could write it here but it's It's a really huge polygon, and I don't want to look at it. You can do it at home. Try this. It's a huge string, so I could actually do quickly len. No, it's not possible because of dictionary. I can't write the len. I could make the polygon, but it doesn't matter. So let's look at name. It's Indonesia. It's just the first entry. It's just completely random. And let's write the name in Chinese. and let's look Indonesia, which continent Asia of course and let's look at the population and in which year was this population recorded so 2017 they had 260,580,739 people there so it's great information so let's go through the whole data set and let's just look at the Germany polygon at the Germany entry and I want to know how many people live in Germany and I want to know what type is this geometry and I want to store I just store in geometry this name this variable geometry I store the polygon so we see yeah 80 million something I guess that's probably true 2017 multi polygon multi polygon by multi polygon yeah there are some islands and yeah exactly exactly so it's correct it seems later we will display it but first we want to do some we know two data sets now we know this airport data set and we have a polygon and we know there are some polygon operations so why not do something very simple i create this multi polygon from i create a shapely object out of this multipolygon object. I call it DE, Germany, doesn't matter. And I take a point and I check if this point is within this polygon. It just writes this airport. So let's look what comes out. Oh, that's nice live coding. I did not execute that. It's too difficult to execute, I guess. Okay, you see we have some airports here. What I have to mention is this airport data set is not complete. It's public domain data set. It's not official. Some people weren't here to do that. So some airports are missing, and we see this list. But it works. It's actually testing the polygon. So you see, the great thing is in about five, six lines of code, You can do such operations. And as I mentioned before, if you programmed this yourself at one point, you would be crazy to do such operations. It's really simple. And we didn't go to GeoPanda yet. You can do the same with raster data. I'm actually skipping this. If you are interested in raster data, you can try it yourself. Here we open some GeoTIFF, and we will display it, rgp channels but i'm obviously running out of time so i'm going to geopandas now so i open a data set it's called cities 5k this one i want to open to show you cities 5k csv you see it's csv and in in jupiter lab it's opening nicely you have all cities of the world which have a population greater than 5000 in this data set it's really it's not a huge data set, but it has more than 50,000 entries, so it's a nice data set, actually, because you have many properties, many metadata, time zones, et cetera. This one comes from GeoNames. It's also open data, and you have, of course, the center coordinate of the city and some things. So, I guess you know pandas. What a silly question. And so we opened that with Pandas. I do separate. There's no header and so on. I do low memory falls just for fun. And I want to see the data and to see if it works great. Somewhere we have the coordinate. That's geographic coordinate again. And what I'm doing now is I'm just reducing this data set because it's too many properties. I'm just interested in the name. latitude, longitude, and the population. You see 0.0, there can be half people. So I could have done it better, it doesn't matter. So we can do things, we can do queries, we can check cars, yeah, it's actually inside. So I don't know if this is true, 799.0, one is missing. So now, you see, we could do things, But the problem is with pandas, the geographic information is just numbers. And there are some things we could do like measuring a distance to something or new queries. And for that, GeoPandas was invented, which is basically pandas plus geo. So what we have to do is we have just to create a column called geometry and add that at the geometry there. And the geometry column is just a shapely point. See the shapely and I just do it quickly Should This connection here is not too stable doesn't matter at the moment and we see we have a geometry column and here is actually the WKT string inside and This data set we just have points so i can remove the lat long i don't need this anymore and i just call data frame plot and you see if i plot all these cities with more than five thousand oh i forgot to these are the things matplotlib inline i have a notebook server where this is pre-converged matplotlib inline so now we see it In an instant and these are all the cities greater than 5,000 you see it It's a nice map of the world more or less with some some white Errors where there are no cities with more than 5,000 But you guess it seems to be correct. I Can export that to a to a common vector format? I'm not doing that and I'm I'm just storing this city for later And what I'm doing now, I create a big cities data frame with population greater than 10 million. And you see I'm sorting this too. And you see this data set, Shanghai is the biggest. As I said, it's an open source data set. It's not the best data set, but you get the idea. If you use a real data set, it works much better. But this is what we get, open data at the moment. And of course, we can load this country data set. and we see here, I do it quickly again, I just take the most important columns, and we see if we take the country data set, we have the polygons in the geometry. What I did here is just loading the shape file, that's the vector format, and I get the Germany polygon out of it, and I just plot it, this GeoDataFrame plot, and we see, if you know the shape of Germany, It seems to be correct. And here you see why it's a multipolygon again. This is these islands here, which are separated from the main polygon. And I can do the rest of the world, all countries which are not Germany, and then plot both, black and green. Don't ask me about the colors, just random. And I can check the geoseries here. I'm going a bit faster ahead. I can, of course, check all cities within the Germany polygon again and create a data set out of it. See, it's running. This was inexpensive. The Germany polygon is too complex, I guess. And you see, oh, I didn't know. La Calmin is inside Germany. It's in Belgium. So this polygon is just at the border, and so it's not perfect. If you have time, which I doubt, we will see that this polygon is not inside. So let's go quickly to another library called Folium. It's based on Leaflet. It's a JavaScript library. And this is a very simple way to create maps. But because now it seems I don't have internet, no internet, it looks bad. I should just, why I have, I should have the network here, you see it seems to work now slowly, let me execute it again maybe, what we would see here if the internet works, ah I should not put this outside, but it seems to be a problem. No internet. So, yeah, I guess here is the normal map of OpenStreetMap. Because we are running out of time, I can't show you all details, but I do have an already executed version which is still cached, so I can show that quickly. It takes a few seconds to open. How much time do I have? Three minutes. Perfect. Including Q&A? No. Okay. I used network care. But it doesn't work, it seems. I mean, I tested before it worked, but it seems the VPN is probably the problem. It's possible. But let's just scroll down. You see here, I executed this cell in the before, so you see here how nice, you see, it works. I just plotted here the polygon from before in different ways. And you remember the big cities, maybe you remember this big cities geopandas frame with the population bigger than 10 million. I can just add this one, create markers out of it, and you see here the markers of the city. This was too quickly. I think now my computer starts crashing because this notebook is too big. And here you can create different icons. You can click on it, and you see the population on it. That's quite cool. You can save this as HTML, and then you have created a map of this with this. The problem of that is you can't create dynamic maps, so moving objects and so on. And for that, there is another library called IPyLeaflet. It's actually a Jupyter extension for the Jupyter Notebook and for JupyterLab. You just install it. Actually, when you want to install it using the lab, you have to install Node.js, and you have to install the JupyterLab extension manually. That's a little bit of a problem. This takes a while, but once you have done it, you can do some dynamic maps. So you can have interact. You see I can't display it, unfortunately, because no internet. But you can imagine, you can change things, we can move things around and so on. You can try it at home. So I get the stop sign. You can take time out of the Q&A if you want. I don't take time from the Q&A. I want to show... No, okay, it's okay. We go to the Q&A. Thank you very much for your attention. Thank you a lot, Martin, for the interesting talk. So let's hear the questions.

Speaker 2 [26:07]

Thanks for your talk. I have a question related to your example where we saw the long load time for all the cities within Germany face a similar problem. And I feel like one could use bounding boxes or other bounding elements there. But is there something implemented in any of the libraries you showed? Like compute a bounding box for my geometry and then check whether it's in there?

Speaker 1 [26:35]

Actually, of course, in this short time, I couldn't show you all the functionality of this library.

Speaker 2 [26:35]

Actually, I...

Speaker 1 [26:41]

But in Fiona, you have this functionality, just bounding box, and you get it. It's no problem. And also in GeoPandas, you can get the bounding box, one command, and you have the bounding box. This is all implemented.

Speaker 2 [26:51]

Also for multi polygons.

Speaker 1 [26:53]

For multi polygons, for multi lines, for multi points, for mixed things, anything. It works.

Speaker 2 [27:00]

works. Okay, thank you.

Speaker 1 [27:03]

Exactly, and then you will see it. Exactly. Any more questions?

Speaker 3 [27:18]

Thank you for a nice talk. So I was working with Geopandas and I was matching like 17,000 points with 500 polygons, which were police districts in the Czech Republic. And I found out that it either doesn't scale well or I did something wrong. What's your experience with scaling?

Speaker 1 [27:38]

Yeah, actually, what I didn't mention, and it's good you mentioned that, the problem of pandas, and let me also say geopandas, is you can't have, it's all in memory. So you can't, it's not chunked in any way at the moment. There are plans to do that. The people, the authors of geopandas are working on that. It's not that easy, but they try. But at the moment, if you have this many polygons, then I would recommend using Fiona directly. And like this example I showed, where you can match if point inside polygon or polygon inside polygon or line inside polygon or whatever, you can do quick two loops and check this. It will be 1,000 times faster. Two more minutes for questions. All right. In that case, let's have another round of applause.