Geospatial Data Processing with Python: A Comprehensive Tutorial

Speaker 1 [00:02]

Thank you for the introduction, and welcome to Geospatial Data Processing with Python, in Python, with Python, a comprehensive tutorial. When I wrote this abstract, I was very optimistic, and as you see, I wrote many modules here. However, I gave this tutorial already one time this year, and there is no way we have time to cover all. So I highlighted the most important ones and in the case we really have time I can Give more information about the other ones, too All notebooks are available Under this link here So you can find this link in the notebook typical chicken egg problem and Those people at home I have a second screen here with this call at least I tried however my disk got just crushed and is downloading some updates at the moment. I do have Discord on my mobile phone so I can see your questions so feel free to ask if you're at home. And the first question already came, will the notebook be shared? Yes, you see it on screen the link and you also see the link in Discord. So let's get started. We need GeoModules and if you install Anaconda you can do the usual stuff. However some people prefer to install things other ways so I provided basically three let me say four variants how to install all GeoModules. The most easy one of course is to download Anaconda. You create your environment and go on. Some people don't like Anaconda because it's big, you need disk space, if you only do your own environments then you waste disk space in the sense you never use the base environment. I myself use Anaconda, at least not on this computer but on my computer at home. Sometimes I prefer Miniconda, it's basically the same, just doesn't come with all the things pre-installed and some people they say Mamba is the best, of course if you want you can use Mamba. Then there is one fraction, I don't ask the audience, who still are convinced pip is the best, so be it, you can try it, there are some binary modules available, especially on Windows, on Mac you will run into trouble, on Linux you may get it, you can do it, it's possible but it's really not worth the hassle, just install MiniConda and you are done. Then to install our geospatial environment you just follow these three instructions. We use everything from CondaForge, now you may ask why you could omit this CondaForge, but some modules here are not on the normal Conda distribution, so I really recommend using CondaForge for all modules. If you do a part without CondaForge, another part with CondaForge, then you get a pretty messed up setup, so follow this. I use the convention that I use the Python version in my environment name, so for example this environment I call GeoPython 3.10. We have higher versions of Python available at the moment, I'm aware of that, we also have lower versions of Python available, sometimes you can go back two, three numbers because some modules are a little bit behind. At the moment everything runs with Python 3.10, so that's fine. After installing all this stuff, it takes a couple of minutes, you are ready, don't forget to always activate your GeoPython environment and then you run JupyterLab and you should see the same thing like I do now, JupyterLab. I know it's a Python conference, I don't really have to explain all these steps, but who knows maybe someone still needs it. If you are seriously using GeoData, I really recommend using GIS software. There are tons of GIS software available. I provided a list here, best GIS software, so I'm not going here and you can, my net is I think my internet is not working properly luckily I prepared okay it's working connection was interrupted that's probably the reason why this cord is failing to update so okay so I recommend using QJS it's open source it's free you can get started with that and it runs on all major platforms anyone using Linux some people Linux no problem it runs on the Renault everything here is tested and runs perfectly on the Linux anyone using max Oh more use max than Linux interesting so no problem too on Mac OS you will not run into problems, you can install QJS, everything works, you can install all packages, and I assume the rest is using this, how is it called, Windows, yeah, okay. So I installed QJS, I will not cover too much of QJS in this tutorial, I will use QJS just to open some files maybe to show you it really works, that's basically the thing I will do now in the next couple minutes so that's the notebook number zero and let's go to the notebook number one if you work with geospatial data the most important thing of course is geospatial data so you need your source if there are two ways basically you can obtain data actually three ways you can obtain your data from an open data source. The second possible solution is you pay for this data. There are tons of really high-priced data sets available around the world. Every country basically sells its data nowadays. Some keep it secret, some sell it. Some are open data. Some are open data, for example I'm from Switzerland, in Switzerland we have everything available as open data, at least base data. We have 10cm orthophotos of the whole country, downloadable, it's a couple terabytes, so we can download 6, 7, 8 terabytes of data, completely free, and use that. We will not download 7 terabytes of data, don't worry. The third option, if you want to have some geospatial data, is you create the data yourself. You can buy a drone, for example, you can fly over some areas, actually nowadays in Germany you need a drone license, so you have to pay 300, 400 euros first to get the drone license, then you are able to fly around and collect your own image data. Or you can also create vector data yourself, for example, using OpenStreetMap. But let's first not do too much with Python. Let's do two things. First download this 309 megabyte data set. I'm laughing because this will be fun. So first thing, I create a folder called GeoData. Luckily I already did that in the hotel. It was not best. I realized yesterday that I don't have the data with me, so I downloaded it myself this morning actually. And there are two functions, download and zip, which is in this download.py file. Why do I do that? Because I love progress bars when downloading, so this function I created a couple years ago just writes a progress bar while downloading, unfortunately I can't show it now because we don't have internet obviously. So you just download it and as a data set we use this Natural Earth data set, it's available here, I can't open the link unfortunately, they provide some free vector data and some raster data, it's a great source to get started with Geo data, however it's not the best data. So for example, they have airports data set with some airports. It's great, but it's not complete. So it's a nice data set to get started with the formats, but I would not recommend this data for everyday usage. There are, of course, other sources of geospatial data. There are some links where you can just Google actually it for your country, and maybe if you are very lucky, you get some open data. So the first thing, we download this 310 megabyte data set here, this raster data, and this function downloads it and this function unzips it. I did not import this file, so it got narrowed, don't worry. So I already did it, so I have a flag here. If this is set to false, it is forced not to download it again if it's already downloaded. So because I already downloaded it, I get this message, file already exists, not overwriting. The same thing with this unzip function. This unzip function only unzips files that are not yet unzipped, just sometimes I really have huge data sets and if I re-execute the whole notebook, I don't want to wait 20 minutes again if you download that and you get an error because now we have a very bad network then you may have to delete the file manually because it's corrupt file so just in case but unfortunately I think did anyone manage to start the download one hour left perfect that's actually why I'm I did put this that's why actually I put this notebook already the link on the program already three days ago if you really seriously just kidding I'm professor I'm a professor I you are not my students okay okay so this will answer bit of course I would recommend doing it at home so I don't think this will work here too Actually, usually I use Blue Marble data set from NASA, which is around 1 gigabyte to download, so I especially chose a smaller data set because I knew the internet will not work properly here, but unfortunately this is what it is, and bringing USB sticks would not really help anything in this group, I guess. So you have to trust me, this will download. see a nice progress bar. You see the progress bar? Slowly, slowly. Okay, great. So I will ask you in an hour again. So once this is downloaded, one thing, we have the GDAL. Unfortunately I really don't have internet, which is really a pity. I was hoping I could go to the GDAL website unfortunately this does not work what I wanted to show you is the GDAL library GDAL is a geospatial data abstraction library who of you already knows this library so many not all not even a third okay this is now a num focus project I heard since a year I was actually not aware of that I just saw it it yesterday in the opening session. It is also an OS Geo project. It's a very old library written in C++, and I used it in the early 2000s the first time, so it's really a very old library. And it's really great because it can read all common raster formats. And what I wanted to show you, if you go to gdal.org, if the web page works, you will see a list of all supported raster files. For example, you all know JPEG, I guess. You probably all know PNG and all these formats. There are more formats, especially for GeoData, where you need some additional information to GeoReference. So you want to know which position is my pixel. And for that, we need more information. We come to that in an instant. But let's go to another thing. GDAL also comes with a ton of utilities. Some of those utilities are written in Python, and some of those utilities are written in C++. So you can just call that. And in this first step, I just want to show you some, just three maybe, of these programs. One is gdal-info. If you don't know this, with this exclamation mark, you can call a program within Jupyter Notebooks. So for example, if you are in Windows, you could write dir, and you would get the list of your files in Linux or Mac OS, you can write ls, and this actually calls the command ls, and you see the output of the current directory, you can provide some options, list all here and see, very secure. I'm a Mac, not on Linux, so that's fine. And you can call commands with this exclamation mark. I think most of you already know that. Who didn't know that? Now I have to ask, who knows that already? Okay, okay. Not all, so it's always good. You can also do some stupid things, by the way, if you don't know that. You could, on the Windows you could do this, for example Calc, you can try it and then the calculator opens. Don't do that of course, it's also a reason why some companies don't like Jupyter Notebook because you can really do stupid things, for example such things, no no don't do that. Let's not go to that. Let's go to GDAL info, I also teach cyber security by the way, so okay, GDAL info, so I provide this TIFF file, TIFF is a known image format, there are some variants in TIFF, there is, you can actually, basically you can, TIFF has a header where you can write your own stuff. So there is a standard called GeoTIFF, where some geospecific information which is standardised in the header. So we see the driver is GeoTIFF, so we know now this file is a GeoTIFF file, and if this gdal.org site would run at the moment, you can look it up, GTIFF, and you see all the options you have for this format, unfortunately I can't show it now. You see the file name, of course, and then you see the size, and now you know why it takes an hour to download. It's 21,600 by 10,800, which is really a small image in the geodomain. I really took a small image. I managed to download it using Hotspot on my mobile phone, so I don't have to do that. Okay. So I heard that's a good... So you can download it on your mobile Hotspot if you have one. I actually don't have 5G here, but if you are lucky, you have a fast, maybe it's fast. It's really, actually yesterday I, yesterday I used the internet here, actually only upstairs, maybe it's downstairs here the problem that, it's a pity. So that's not good. So that's the feedback we will give to the conference. Okay. So it really does, you can look it up at home, this $G site, I really recommend it, you will see the list it's really easy to see and now okay so now it wants to password show password it's very safe oh now see it that's actually it's recorded so that's a other say cybersecurity so the next thing you see here is the geographic coordinate reference system and this is wgs 84 wgs stands for world geodetic system 84 stands for 1984 that's when this system was created this is used in gps so we have geographic coordinates and we have actually a ellipsoid which is quite which is quite good approximation of the world as you may know the world is not flat so we need an ellipsoid to approximate it and you see here some parameters of the ellipsoid, this one is the radius and this one is actually the second radius hidden, it's flattening so the earth is a little bit flattened towards the pole, not flat, flattened, not the same. And we have also length unit, this is actually wrong, length unit is in degree here, and we can't cover that. So we see origin is minus 180 degree longitude, 90 degree latitude, so the north pole would be here actually, that's a point and here we see the pixel size in x direction and the pixel size in y direction in map unit so it would be in degrees of one pixel is 0.016666 degrees that's very important to know with this information we can really transform pixels to geocoordinates so we can open it i will not do it because i provided a screenshot here we can open this tiff in qjs just drag and drop it here and you will see this wonderful world map which we downloaded it from natural earth why do i open it with um with qjs not in the notebook yet because it's tiff and tiff is usually not supported by web browsers web browsers can open pngs jpegs GIF and some other formats maybe, but usually not TIFFs, so you have to look for another way to display it in Jupyter. There are also other ways to georeference such images, for example using world files, so sometimes you will see TFW files, actually if you download this and it's unzipped, this natural earth it's here. There is also a tfw file and in this file you find this information it's plain ASCII. The first line you may remember this value that's the pixel size in x-direction, pixel size in y-direction and here is the pixel position of the center of the top left pixel. So with this information we can do a fine transformation and we have every we can transform every pixel to real coordinates there is also a sometimes a file called P or J which is the coordinate reference system so we know it's in geographic coordinates this is a little bit outdated and should actually not be used used anymore, but it's still used. There are more modern approaches to do that. For example, use GeoTIFFs and use the image header or use these XML files which provide all these things. You can use the tool GDAL translate actually to convert outdated file which does not have this information using GDAL translate and you would say, this is not necessary here because it's GeoTIFF already, but you would say, I take my, for example, JPEG, I translate it, I want the output format, for example, GeoTIFF, I want to compress it using JPEG or deflate or whatever you like, then you would say, okay, my spatial reference system is WGS84. That's written with the EPSG code 4326. I will come to that at a certain point probably. And here the bounding box. So that's basically how we create these things. There are other applications I provided some examples here. I'm not executing that. You can do that. It takes a couple seconds. So another thing is we can change the image size using GDAL translate just by providing the flag outsize and then you provide the new width and the height you can also set 1 to 0 if 1 is 0 then it takes the correct aspect ratio that's sometimes a good thing because maybe sometimes we are too lazy to calculate it ourselves so we don't know the aspect ratio so this will provide the correct aspect ratio you can also do with a percentage, so we can say it's 50%. You could also set the second one to 0, then the aspect ratio is taken again. So we can try this out quickly. You see, no, that's a problem of Jupyter. You see, this takes a couple, less than a second maybe. And we have a new small version. And we see the new small version is still a GeoTIFF and is now 2048 by 2024. Great. So we can also extract, and that's something which is actually quite good. We can extract a rectangle out of the data using GDAL translate by specifying a projwin. Projwin means I used the coordinate system of the image, so in our case, it's degree. And I'm from Switzerland, so I take the bounding box Switzerland and store it here and then actually this this new TIFF file is let me open that in QGIS is actually the extract why is this important as I said before sometimes we do have no that's wrong actually I have to open it here great now QGIS is somewhere hidden I don't have a mouse and I'm not good with these trackpads as you may see so this new file Switzerland TIFF it's 1055 so you see it's now it's not some magical trick and we see here the extract in this is a Switzerland and I won't do the image is small so why is this important sometimes we do really have huge data sets in the terabyte range and you just want some part of it to actually open it. You can't open that easily in QGIS. Actually you can, but I'm not going into that today. Easily in QGIS, so sometimes some extents are important. Now we covered RGB data, so we have red, green, blue. Color image is raster data, can be more than that. There are, for example, many use cases where you have satellite images which are multi-spectral, so you have maybe 12 bands from infrared to I don't know what. And there is one other thing, it's elevation data where you actually have one band, not RGB, just for example floating points or integers where you have elevation values in that. And for that we can download, unfortunately it really doesn't work, this link. Ah, it works now. Okay, great, so maybe I can also open Discord now. I'm not GitHub Discord, these logos are quite, okay. So here is the data from the shuttle mission from NASA SRTM, you can download some tiles, you could click on all if you want. You can select some, I actually did one time click on all and download it all, so I don't have to do this anymore. You can search and you get actually some, you see a preview here and elevation data you can download here. I clicked one time on all and created a small script. Because this download links, I inspected the page source and then I extracted these files. So I downloaded all. The system is quite easy, it's the file name here, you can do that. I have a backup server in case the first one doesn't work. So you can download and unzip it, I already did that. And if you really want to download everything, I provide this link here, don't do it now, it's over 30 GB, so you have all the elevation data around the world. But please don't do it now, that would be even more congested. So once we have the data, now Discord is actually working again, so I can answer the questions at home if there are any. No, there are no questions, I see you provided the link, perfect, so now you can, also here you can ask questions anytime, just stop me, it's meant as a tutorial. So the next step is we can create something called a hillshade using the GDAL DEM tool. And this hillshade is just, if you take elevation data, it's just a nice view. You may know this from some maps you already saw. And let me, what was the file called, srtm underscore hillshade. So it's this one. And you already see it here. I also took Switzerland because of the mountains. I didn't want to take Berlin because there are not too many mountains. So you see here, you see two things here. The first thing, it actually fits. The other, see, its position is correct. Because all these files are georeferenced in the geotef header, what we saw before. And we see some mountains are shaded. So sometimes this is nice. There is also a second variant called color relief map, where you can provide the first value, you create a text file, the first value is the elevation, so from 0 to 500 we have this RGB color, from 500 to 1500 we have this RGB color, you could do a really big file. I create this text file in Python quickly, so this is a normal text file, color txt here, color txt here, see, normal text file, just created in Python, because of the fun for it. I like notebooks where you can just execute everything and everything works without creating something yourself, so I create this text file, which is completely stupid, but I did it like this and then we make this color relief map you see this takes a little bit more time it goes through the whole file what we didn't see is with GDAL info this file is 6,000 by 6,000 pixels usually it's faster but I think my notebook is not charging so I need an emergency help from the I think I need my power adapter this power adapter does not work so there are actually more cheetah programs while I'm looking for the power I can show you this call relief map here in QJS and you see the wonderful map this is actually yeah it's slow because I'm running on battery you see this yeah Bands? No, bands are red, green, and blue in this case. In geo raster images in the geo domain, we speak of bands because RGB is just a special case. I mean, most or many data sets are not RGB. They are multispectral, so you can have 15, 20, 30 bands, and this is not RGB. We don't speak of color or something there. And for elevation, we usually have one band with the elevation value. So, I need power, I think this is long enough. Okay, perfect. So, the problem is probably solved. Another question? Yeah, question. So, all this can be stored in a TIF file. Like, if it doesn't have some kind of limitation, like the bits or how it's... Nowadays, TIF doesn't have a limitation anymore. If you have older software, you you may run into problems with TIFF that are bigger than 4GB. There is a standard called BigTIFF, which is included in TIFF nowadays. So TIFFs can be unlimited, basically. Is SRTM data the only one that has the elevation band? You should use the microphone for the audience. Check. Is SRTM data the only one that has the elevation band? Or is that available in other packages? No, there are actually many data sets around for elevation. Most countries, for example, have the data sets. Some are free, some are not free. I'm not sure if Berlin has a very nice open data policy. You could probably download Berlin elevation. It's not interesting. There's not much information there. But it's just very downloaded. TIFF is one of many. Now I can actually show you the GDAL web page. maybe not, with all the raster formats. Actually, this is the wrong page. With all the raster formats, if you scroll down here, you see all the raster drivers. And you see this list here. This is all different raster data sets. So for example, we have JPEG 2000, which is a wavelet-based format. Or we have some completely unknown formats. We can also have PNG with additional files, so it doesn't really matter. However, GeoTIFF or TIFF, in TIFF, you can really store multiple bands, not all formats. For example, in PNG, you can't store elevation data in float 32-bit. So for that, you need a format which is supported. But it's not only TIFF. There are other formats, too. OK, we are running out of time, as I feared. Let's go to the next thing. GDAL comes with, let's go to Python now, finally. GDAL comes with a binding, Python binding. However, GDAL is written in C and C++. And it's this, if you know C, it's those structs. And you have the file format as a struct. So for example, float is three, int is two. And if you do a binding in Python, it's horrible. So you basically have a C binding, a direct C binding from Python, and you can work with it, but it's a mess. For example, in Python, we are used to this file. We have some file, open. So the first thing is we don't use a capital letter in function names. In GDAL binding, we have to open with a big O, and it's just not nice. So some people created a wrapper around GDAL, which is called Rasterium. And Rasteria also has some additional things there. I will not go into details here, but Rasteria is just a better way to use GDAL. But Rasteria uses GDAL, you know, the principle. So it's just a nicer, sorry, I'm losing my voice again. So we import Rasteria, we import some matplotlibs, we want to display these things later in matplotlib. So we can just open our dataset, this is the smaller one, I hope I executed this one, yeah. And we can look at some metadata here, and that's a thing that is stored in the header is now displayed here. We have a geotiff, we have a data type, which is uint8. This is actually from NumPy. So we have 8 bit per pixel, which is normal for images. And there is something called no data value. Maybe some pixels are not covered or some points are not covered, so we can have a no data value. For example, say 0, 0, 0 in this image is no data. It's not a valid data. Or for elevation, this is usually set to a minus 32,000 or something like that. This is not set here, however. And then you have the width height. We have the count here, and that's the number of bands. So in our case, we have three bands, RGB. And we have the coordinate reference system. This, again, is WGS84, so we have geographic coordinates. And we have the affine transformation here. Again, this is the pixel size in x direction, pixel size in y direction. And here at the top left, but not the center, it's a top left coordinate, but really the top left, not the center of the pixel. So it's exactly minus 180 and 90. And you can also access this data here. I'm not going to do that again. Up to here. And there is something interesting about this transform. We see here this is an affine transformation, and you can also have the inverse of that by using this tilt operator. And you see this is the inverse, so this one would be the center of the, this one is the coordinate in pixels actually here, and we have just the inverse of this one. We can also display the bounding box of this. We have the full planet, as we see here. And now this is affine transformation. We can transform 0, 0. The pixel coordinate 0, 0, that's the top left pixel. And we get the geocoordinate in degree, in our case, minus 180, 90. We can do the inverse. We can take a geocoordinate and get pixels. We can use another geocoordinate, for example. That's where I work. And we see that's a pixel. We also see, you may say, pixels are integers. Yes, you are correct. However, this would be the correct pixel in between. And we would have to create integers. And then we can access this pixel. So let's store something and get this pixel. We will use that later again. can remember px, py is this pixel coordinate I transformed. So we will see where it is. And now, I actually did that already above, but let's do it again. Import matplotlib, we want to display our raster. And that's quite easy. This would basically be the whole program again. I open it. I used a small one. Actually, I think I didn't execute that before. Sorry about that. I skipped this, never skip something. I skipped something, but it doesn't matter. We can use the big one here. This world tiff. World tiff. Little bit, real time coding doesn't, it's not bad. This will read all the bands. In our case, we have three bands and we read RGB and that's also important. These three bands are usually stored as bands in the TIFF format. So we have really the red color, the blue color, and the green color separate. And we can stack it together using DStack from NumPy. And then we get actually RGB triples for each pixel. So let me show it. R would be just red values, and the stacked version RGB would be just these triples RGB, so quite easy. Most people don't know DStack, but it's exactly for that. So we can use that for matplotlib, for example, we can use the imshow and provide an interpolation. Let me use the nearest. is actually wrong, it should be much faster. Maybe because I took the whole image, which is quite big, not a small one. And we see the point I calculated before, that's where I usually, where I live. Just for fun, you can create another coordinate if you want. There is also a way to use directory as to your plot data, you would just call rasterior plot show, I'm not executing it because it takes a little bit of time, And at the end, don't forget to close the data set again, but this one I will call. So now the data set is closed. However, RGB is still set, of course. RGB, not RGB, RGB is still set, so sometimes it's a good idea to manually delete these things, especially if you have big data sets and hundreds of notebooks, and at a certain point in time, you will run out of memory. It's no problem here, but it would be a problem for larger images. Okay, let's go to the next topic. We don't only have raster data, we also have vector data. Vector data is a very important thing because it's more common than raster data. You can really download tons of vector data sets, which are publicly available. There are a couple of things you have to know about vector data. There are some certain kinds. There's a standard code, OGC simple feature access, where such things are defined like polygons. What is a polygon? I can't draw pictures here, but you can imagine if this point is overlapping, that's a non-simple polygon that's forbidden, for example. can't have overlapping polygons and these things and this is all defined in the standard if you have overlapping polygons you have to create multiple polygons it's called multi-polygon actually and contains more than one polygon another specialty is polygons with holes this exists for example many countries have this with their holes or multi polygons with multiple holes and this is all defined you can create these things and there is one very nice library called shapely which can handle all these things I don't have time to go in depth shape it but for example you can calculate error of such things very easily by just calling error and get error of such polygons even with holes and everything which is really great thing you can also check if there is a point inside a polygon if there these tests and it's a really great library it's based on the geos library geos is geometry engine open source geos it's also written in in C++ and there are bindings for Python available. And also bindings for, yeah. Excuse me. Once I was calculating area. Could you please use the mic? People at home don't hear you. Once I calculated, and thanks for this kind of thing, but I had a file and then I calculated area with cookies and then I compared the area with Python and there was a fraction difference between, I checked the projection system was same. I was making sure that I've repaired the geometries of both. But there is always this fraction which I cannot explain and understand. I can't tell you, QJS also uses GHS. So it's the same library, so this should not happen. So my guess is that something happens in your Python code, but you can show it to me afterwards and I can look at it. I don't have it with me, but maybe I can write you later about that. Exactly, and I can look at it. But it should be the same because it's basically the same library. And that's also a good point because the numerical stability of such things is very difficult. If you have very small polygons, for example, or you have a collection of really huge polygons, and it also depends on the projection system. if for example you have really large numbers in your projection system there can be numerical instability with these things I programmed all this stuff myself many years ago and I know it's really really hard to get the numeric stability for these things and it can be a mess. Which one would you trust most? I would actually use, I would go with geos because everyone else is going with geos too and you have the same results some some data some countries for example country borders and this thing or eras they store it with the data along so you will not run into trouble there so the era of a country is exactly what it is in this data set and if you calculate it yourself you get another result that's quite normal and And there was another question on Slido. Is there any useful Python lib to generate a directory-based vector, also known as XYZ tiles, like tip canoe? Actually, the last one I didn't know. A vector format as XYZ format that would be more a point cloud and not a vector data set. You have to provide me a little bit more information than I can. Who asked that question? Maybe we can go more into detail. Maybe someone at home. You can also use the chat, by the way, in Discord. I have it open here on my second screen. That's quite handy, a second screen for these things. OK. We have these things, point, line, string, polygon, multi-point, multi-line, string, multi-polygon. Then we have the geometry collection, which can be anything combined from these types. In Shapr we do have more things. We also have a box, for example, just for convenience, so we don't have to create a polygon every time. So there are many popular file formats. If you download this data, you will probably run into something called shapefiles. It's a very old format developed by Esri. It was a long, closed format. It is a little bit more open now, but you should not use this format anymore today. If possible, there are really some limitations. It's really an old format, and there limitations which make no sense today. For example the field length is really restricted and only capital letters or you will run into problems if you have UTF-8 sometimes and these things which and at the end shapefiles are a collection of many files so it's outdated format there are better formats also if you if you are using ESRI software you can use GDB from ESRI too and that's a better way to go. Then there's GeoJSON, there's KML, known from Google Maps, there are some other, it's just a collection, you can find more formats actually here on the GDAL web page. If you scroll down you see the vector drivers and all these vector formats are supported by GDAL, which is amazing, and my favorite format nowadays are these two. If you are an ESRI user, I would recommend using GDB, and if you are more into the open source, if you're more in the open source community, this OGC GeoPackage is really a nice format. It's an open format, and it uses database, actually, SQLite databases behind it, and it is standardized, and it's really great because you can have many things in this package, and it's really fast. So let's download some data again, now the hour is over now, not yet, but now it's time to download some other things, luckily I already downloaded it. This is again a couple hundred MBs, it's not big but it's probably too big for the conference again. unzips it here and you will get this folder packages and this GPKG file and here we have some sample data. Again if you are interested you can go to the GeoPortal of your country, you can Google data, you will probably find it and if you are lucky your country even provides it for free or your city or wherever you obtain your geodata and as i mentioned before berlin is perfect for geodata they opened up quite a while ago they even have a very nice 3d model of the city which is freely downloadable which is quite great and inspired many people and scientists for For new projects, I always promote open data. For example, you can calculate the solar potential with this 3D model, and some people can do new things, and open data gives us tools to create new applications. So there is one very nice library called Fiona, which uses GDAL again for vector drivers, And you can see it doesn't really use all of GDAL, but it uses these formats at the moment to load vector data, and the most important, of course, unfortunately, shapefile, and geopackage is in here, and some other formats which are there, and we can, of course, test if something is supported, and yes, geopackage is supported in this version of FIONA. So let's list our so-called layers from this GeoPackage, and we see this is the natural Earth data set. We have many layers, we have boundaries from disputed eras, unfortunately this exists in our world, and we have countries, we have railroads, we have lakes, we have really nice collection of vector data here, and we can try to read it, it's just calling Fiona open, specify this layer. Actually layer are something specific of the geopackage. If you use shapefile, one shapefile has only one layer so you don't have to do that. You would have for each layer here your own file. And the next thing you probably want to do if you don't know the schema of this file, you can see what is actually in our airport layer. And you see we have properties and all these properties are in there we have a scale rank we have some type what it is we have names in different languages there that's called properties every vector data set has different properties of course and the second thing the most important thing basically is the geometry and in our case our airport we do have points so for every airport in this data set we have a point and for every airport in this data set we have these properties and we can look at them. So we can, of course we could just get them separately over typically Python, it's a dictionary, you can call the properties here. And of course we couldn't do the same geometry and you see it's a point. In other layers it may be a polygon, it may be a multipolygon or something that's specific to our layer. And then we could get the first entry by just iterating through our C, C was remember we opened that, this Fiona open, called it C. And we can iterate through these airports and I call this variable airport and I get the property's name of this one. And we see the first airport in our data set is this Indian airport. let me not pronounce it, I don't know how it is, and we get coordinates. And of course all other properties we want. We could just call all these things, for example maybe the name in Chinese. We could try if this exists. Okay, it's an Indian airport. I didn't call close yet. We could do that. So airport properties And then name, Chinese name, let me see if this works, I don't know actually, that's a fun thing about live coding, never know the result, actually it didn't work. So it's not set, so we have to get another language. Maybe, but it's he, Indian, maybe, so it's Hindi, and we see the set, which makes perfectly sense. Wouldn't it have to store the projection system? The projection system is stored, we will come to that. The projection system is stored in the file. We can, I didn't close it yet, we can look it up. Actually, it should be before, no, we didn't see in the scheme, we don't see it. We actually have to manually call it with CRS. I didn't close it yet, luckily, so I can do that. See CRS, and we see this is EPSG4326, which would be WGS84, so we have our geographic coordinate system like we had before. So it's covering, we have coordinates in degrees. I didn't mention that that's in degrees. Latitude, longitude, we don't really know which is latitude, which is longitude yet, because some people prefer to have longitude first, other than latitude. So that's sometimes a problem. Let me speed up. You see here I can just print the first 15 names of the airports by iterating normally to this using a for loop and I can for example search for a specific airport and it doesn't find this one what did I do wrong I was thinking it's BER Berlin so maybe it's not yet in the data set so let me do another. This one exists. What was the Berlin Airport before? It was BER2, I think. Tegel. Tegel. TXL. TXL. This one. Maybe this. Yeah, this one is set. Okay. This doesn't exist anymore, but it is. Actually, I have a GitHub repository. Someone could report So, we can actually do the same for the countries, we look at the scheme and we see this one is for the countries and the new thing is our geometry here is not a point anymore, the geometry is a polygon of course if you have countries and you could also get the first one which is completely random, the first one is Indonesia. Here we see the Chinese name works, we see the continent, we see the population, that's important and when the population was counted of this country. So we can do that, we can search for Germany and this is Germany. Isn't it nice? So I think it would be great to actually plot this one and there are many things how to to plot this data. However, let's do a shortcut and go to the wonderful module called GeoPandas. You've probably heard of pandas. You may have heard of GeoPandas. GeoPandas is, is Joris here? No, he's not here, so I can do the simple way. I can say GeoPandas is basically pandas with a column called geometry and which stores geometry. That's a simple way to explain it, but it's not that simple, course. It can do lots more things. So I import GeoPandas and now we can just use GeoPandas. I used a short version like Pandas. I think most of you use PD for Pandas because we are all lazy and don't want to type Pandas every time. Same as GeoPandas GPD and we use read file and specified file name and in our case we can also specify the layer let me take our well-known data set with the airports again and now we can actually see our data frame geo data frame with all our nice columns here and the last column is called geometry and we have points here airport and if we open the country would have polygons and multi polygons here and here we actually see more easy it's much nicer we can actually see which which things are set to see some some airports don't have it's not a number here it's not set probably and some some are set and previously we tried this and it didn't work so I usually when I have to beat notebooks I just take the most important columns so we actually see better see better we actually see the geometry here and you have the most important things you see there is scale rank etc and of course we can do all the things we can do with pandas we can have the head we can sort values by column name ascending true false this is standard things we can do queries like we can do in in So that's it, and we can do much more, we can actually plot with it, and that's what we will do soon. But before we do that, we have to cover one more topic, it's coordinate transformations. I said a little bit before, WGS84 and these things, it's not that easy. Usually every country has its own coordinate reference system, some countries share some coordinate reference systems, some countries have multiple coordinate systems. There is a group called EPSG, it's the European Petroleum Survey Group, they collected all, They have an interest in mapping, obviously. They want to know where they find their oil or whatever they are looking for. And they collect all coordinate reference systems and give them a number. So previously, we used this 4326. The EPSG code 4326 is the WGS84 ellipsoid. So it's a global geographic coordinate system, and that's the one you use for GPS. And there are tons of others. I've provided some interesting ones, not all. We have UTM zones around the whole globe. That's, for some countries, very important. And they use those UTM zones as projection system. And you know this one, the 3857 is the web marketer. You probably know this from Google Maps, OpenStreetMap, or whatever favorite mapping software you use. This is the projection system on Google Maps. You know this. You know Antarctica is the biggest continent if you know Google Maps. And we will cover that in an instant. And there is a module called Pyproj, which is a binding to approach 4, again, the C++ library, which can convert those coordinate systems. For example, we can say, I would like to convert from WGS84 to our web Mercator Google Maps projection and get this information here. I do it quickly for this 2056. I'm from Switzerland, that's a Swiss projection system, LV95 it's called, and of course I have to import pipe roach first, and we see from this geographic coordinate, I would get the country coordinate, which is a metric system actually, so this is in meters, so if I go, for example, more north, more east, or whatever, then it decreases, and the unit is a meter, so we can do better calculations here, for example, if we want to calculate distances it's a mess in a geographic coordinate system so it's always good to have a local coordinate system to do that we can also transform back however if you compare these values we'll see we have a little bit precision loss because there is a always a numeric instability between such systems one thing we always always want to do is we want to know the shortest paths around the globe. As I mentioned earlier, I hope I didn't surprise you, the Earth is not flat. So calculating distances is not that trivial because we are on a curved shape. And believe me, I did it. It's really a mess. There is no direct solution to do that. You have to approximate it. And luckily in Pyproch we have a very easy way To do that, we can call the NPTS function, and we provide a start longitude, start latitude, end longitude, end latitude, and the number of points in between we want, and we get actually those points. So let me do that from Berlin to New York, classical route. So we would get, this is the BCC in Berlin, and this is New York, and this is the path from Berlin to New York, the shortest path. if you want to walk there, take this one, or swim, and we can look it up in GeoPandas, and here is a way how to create a manually, a geometry column in GeoPandas, you can use this point from Shapely and zip those together and you get your point, and we have a GeoDataFrame. What I didn't show before, we can actually call plot in GeoDataFrame So let's do that. We plot our line, and that's the way. Here is Berlin. Here is New York. Isn't this nice? Wonderful. You don't believe me? Yeah, you believe me, I guess. What would be nice, of course, would be nice to have our real map behind. And, of course, we can do that. We can just load the layer admin0 countries from this geopackage we downloaded before, and we can open it with geopandas, for example, and we can just plot it. So let me do that quickly here. If I plot this one, plot, you see we have our countries, isn't that great? And now we can combine this one and this one by just storing, actually this uses, you see see it here hidden little bit, this is a matplotlib behind, so we can store it as access object, this countries for example, and then we plot it over the previously access object with a marker size 40 our way, I do some face colors, edge colors and so on, and we see that here is Berlin and here is New York. So now you believe it, probably. You also see this was not a good idea, a little bit big, should make maybe 69, it's probably better. We see here Antarctica, you don't know this usually, I mean if you look at Google Maps you see it differently and we can of course drop Antarctica quickly and then we can, let me drop Antarctica first, I just remove Antarctica first and then I change the coordinate reference system to the Web Mercator. This one is known from Google Maps, as I said before. And then I plot it again. And you see, this is the world you know from Google Maps. This one is a different projection. And this is how Google Maps presents us. Now you can ask, why did I remove Antarctica? Because this is a Mercator projection. And towards 90 degree in latitude, The Mercator projections goes to infinity. So I can show you that quickly on Google Maps. Let me open Maps, Google, com. Correct, you see I didn't even save that. I usually don't use that. Data privacy reasons or whatever. So yeah, here it is. Isn't that wonderful? Antarctica. Oh, it doesn't allow zooming anymore. You can't zoom more out because people were confused because they thought Antarctica is the biggest continent. You see, you can fit all continents into Antarctica, basically, which is completely wrong. There is a better map to demonstrate that, OpenStreetMap.org. Why is it better? Because you can fully zoom out. and you see the same principle you see here, it repeats because I have the same principle and you see Antarctica is really big, that's the WebMarketer projection you also see Greenland is really big so towards 90 degrees it is infinity and this is not it stops here at a certain position this is not 90 degrees, it stops before If you look at Google Earth, you will actually notice that too. I'm not opening it for time reasons. You will see that this is about 85 degrees, a little bit more. And what is bigger is also a hack in Google Earth because it's filled by some triangles there. So this... Now you know why I removed Antarctica, but you don't know how to do it with Antarctica. I could try to not, actually I could show you what happens if I don't remove Antarctica first. This happens. So the trick Google Maps does is, let me transform first WGS84 to this Web Mercator and let's transform minus 180, so this side, and 0, 0. And we get this value here. That's 2 million, 20 million. Actually, you know, if you walk, maybe someone of you already walked around the Earth. How far is it? If you walk around the Earth, which route would it take? For example, you take the equator, and you walk and swim, of course, around. So how far is it? Anyone knows? Let's ask ChatGPT, then. Forty-something thousand kilometers, yeah. No, not eight. It's around 40-something-ish. I don't know the number, but we can actually read it here. minus 20 million plus 20 million if you take the distance between this you have the exact value in meters of course not kilometers then you have to divide by 1000 and then you have this value so see if you do this you have around 40,000 with this WGS84 ellipsoid so what I do now is I just copy paste these values here as a bounding box from minus this 20 million-ish value, minus 20-ish million, plus, plus. And then I get actually this value. This is minus 85. And what I do now, I do the map in a square. And if you look at Google Maps, you always see the map is a square, perfect square. So if we use these values, if our biggest value is 85 degrees, minus 85 degrees and plus 85 degrees in latitude, then we have a perfect square in this macro data projection and we probably don't have this problem that we have infinity in it. We can try to reproduce it in shape. We have this box and we can do some very nice thing in Geopandas. It's my favorite function. It's called clip. So we can clip this rectangle here by 85, minus 85 and we get a new geodata frame called I call it world 2. I can plot it and you see, no you don't really see it but it's clipped by 85 degrees so a little bit is missing and if we convert it now to the web marketer projection we see this actually works and this is how you know the world if you only use Google Maps. So that's how Google does it and the north and south pole is not visible on Google Maps so you can't look it up or something. So if you look for ice bears at the pole you can't find it. Okay, that's it. We have 20 minutes left. That's perfect for two more things. Now we plotted with Matplotlib and Geopandas, etc. There are two very important Python modules to create maps in HTML so maybe you know the leaflet library, it's a mapping library it's actually what is used in OpenStreetMap so if you go to OpenStreetMap.org you see this mapping software which is open source and this is probably done in with leaflet. Sometimes I change, so I say probably done with leaflet. There are multiple mapping libraries in JavaScript available, and volume is basically, it's not a wrapper, but it creates JavaScript and HTML code for us, and we can just look at it in Jupyter, so we import volume, and this line here would create our map, and we see wonderful two lines of code. actually one line of code I could directly write this so in one line of code we have a map with Berlin in the center why is it in the center because here is the coordinate of the PCC and we can zoom in we also so we also set this zoom start to 8 I could increase that for example 12 and then of course it's It starts to increase, so we can create a map in one line of code in Jupyter Notebooks. We can store this map as HTML, can just save this one M to HTML, and then it stores it as HTML. We have some different map types. For example, we can use this watercolor. Then you have this watercolor map. Berlin I guess in watercolor so that's nice isn't it and some other layers there are some some services around I have to admit I don't know the license of this one how if you are able to use that from Esri or if you have to buy a license for that, but it's open, it works, so I don't know, so I will ask one time, we have this nice map from ESRI, and there are some world imagery satellite imagery here, you can zoom in In Switzerland you can use it, this is also available, but the problem of Switzerland is that this coordinate, I don't know why I did this coordinate in here, it doesn't work because Berlin is not in Switzerland. I think I was sleeping while I created this. So I copy pasted all the VCC things and the same is here. So we can add some, actually we can create some layer control where you could actually put all inside, you see here you get this nice icon and then you can actually switch the map of course Switzerland is here so we can't really do that and you see we can switch, we can just add it, some features, I'm not going into volume now, there are tons of features, you can even overlay your own image for example there, if it's georeferenced, you can overlay that's what we look at now we can overlay some things called markers here for example I put two markers for example two hotels we can also customize these icons with the font awesome for example here a beer or hotels whatever and we can store it as HTML so if you create this map we add some markers here We can just store it as HTML and then you can't believe me. I'm not opening it. This is now Also available as a HTML file. You can open it new favorite browser and it looks exactly like that Okay, so now the fun part we can of course draw things here we can for example create polygon here and We see here polygon Just a random polygon we can add pop-ups to polygon by just adding add child pop-up. This pop-up can have HTML in it, so if we click on it, we will see this HTML comes up, hello world, and even image tag inside this, a little bit hidden, I made a normal image tag, and we can have pop-up, so if you click on this, of course you can have multiple polygons, you can have multiple lines, you can have multiple objects here, and then each time you click it, you get some more information, some more information, so you could make some really nice maps with this knowledge. So here, again, polygon, circle, and line. There is a circle marker, there is polygon we saw before, and then the polyline, which just covers all these objects, so we could, yeah, one exercise I often did is with nuclear plants and circles there and see how dangerous they are whatever okay so let's go back to our genetic line from before of course if you have this line we can use I mean maybe we plotted the line here of course we can also try to draw a line or points or markers with with this genetic line and the problem however, is that all web maps, even Google Maps, OpenStreetMap, et cetera, they always have the latitude first and then the longitude. While GeoPanda and many other libraries, they have first longitude and then latitude. And I usually do that with a list comprehension. I just switch this around. If you know an easier way, tell me. I think this is the easiest way actually to swap them. Maybe there's an easier way to do that, I don't know. So I swapped the coordinates here, and then let's do Berlin as a marker, New York as a marker, and then a polyline between the locations, and that's our result. That would be our line. We could also add a marker to the line, we didn't do that. And here, yeah, that's Berlin, and that's, oh yeah, that's New York, okay, so great, isn't it? Some people don't like Folium. There is a second, thank you, I see it, but I'm perfectly in time. And there's a second library called ipyleaflet. And you can basically do the same with it. So I'm not covering everything now. It's just a second module. One good thing about, actually this one is, this doesn't belong here. IPyLeaflet import basemaps. I don't know, maybe it's stored like that. So here, I will update it in GitHub, it's basemaps. And one thing IPyLeaflet has, it has, it comes with these basemaps, and you see these are open or freely available maps. There are some hacks where you can also use Google Maps, but you need an API key officially, so if you want to use Google Maps, you need to buy an API key for that, so that's the official way. Actually, I can't say it, it's recorded, but anyways, it was lucky on a Friday, years ago I programmed my own virtual globe and I had the glorious idea to steal all the tiles from Google Maps. I did that but I downloaded it a little bit too fast and our whole university was banned from Google for 72 hours. Luckily it was Friday and no one noticed that. So we can look at these space maps, it's a, and we can just take one of those, we can just choose the, that's the only benefit of iPilot, there are more, there are some, it's a little bit more interactive and a little bit more integrated in Jupyter to iPilot. I don't really like Jupyter notebooks too much, I mean, they're great, but I don't want to limit myself to Jupyter notebooks. That's why I'm not using these functions with interactive stuff in Jupyter notebooks because at a certain point I want to do my things in a normal Py file for example and that's the reason I'm more on the volume group but at the end it doesn't really matter. You have to know what you want to do and take the best tools for that. Okay, so So we have eight minutes and I have one last slide which I dramatically shortened. OpenStreetMap is a great source for data. And there is one module called OMSNX which downloads street networks and you can also download specific data from OpenStreetMap with it. I'm a little bit covering now the street network. Every city has its street network. I will not use Berlin today. I could use Berlin, no problem, but Berlin is a really big city. And if everyone here tries to download data from OpenStreetMap, we will be banned too. Because there is some usage policy, you can't really download too much from. So it's no problem. If you do it at home, it's no problem. But if everyone here downloads with the same network, then we will most likely get banned. There is one thing I use, so for me it will work for sure, because I use the settings Use Cache. And in the hotel before, I downloaded this already, so it will work for me for sure. But I recommend not doing it here. You can try. If it's too late, it's too late. But I use Basel, Switzerland. It's my hometown. It's a pretty small town compared to Berlin, 150,000, 200,000 population. And because it's small, it's not really fast to download. You see it's fast because it's cached. But if you try it, did someone try it? How long did it take? About a second. Because it's a small city, we can plot this graph, and we see it's really a small city. You can try Berlin. Did anyone try Berlin? You can try it, it just takes longer, but it works, so you have more nodes, more streets. And one thing we can do, we can store this graph, this road network, all these nodes have of course names, there are street names, et cetera, stored there. And we can store it as geopackage, for example, you can also store it as shapefile, but you will get a warning that this function is deprecated and will be removed at a certain point because shapefiles, as I said before if possible don't use shapefiles and we see in this city we have 4651 streets if someone tried it with Berlin this number of course will be much larger. Now something I didn't really tell you before is I used the drive as network type So for cars, there are also street networks for bikes, for walking. You can also download all, and all is private access streets, private access only streets, and you can choose it here. So if you want a bike route, you could just write bike here, and then you can calculate some things with this. Now you are not limited to anything. You could, for example, do the same I do now, calculating routes. That's always a good thing. But before we go there, let me show you that this is GeoPanda's data frame automatically. So we can also just take some street types out of it, and we see we have 2436 residential streets and so on. Actually every city of course has this a little bit different and can do some statistics. You could use Seabourn or whatever plotting library you want to use and maybe display that. That's not that important. You can also export that directly as GeoJSON and put it on a volume or IPy leaflet map. Which is quite funny, I took the CartoDB BlackSync DarkMatter it's called because this is OpenStreetMap data on OpenStreetMap so it perfectly fits so which is not the best idea, but now we can we can call one function, it's called Geocode, that's actually a really great thing where you can actually enter something Basel has two train stations, it's a big confusion actually we have three train stations Basically, we have a German train station, a French train station, and a Swiss train station. But the French and the Swiss one are located at the same place, so it's basically two. And sometimes people confuse these and go to the wrong. So if I go to Berlin, I have to go here. Actually, I can also go here, but that's another story. so I can geocode it and what I get is I get this position back we can do that for basically every every address we want or every location we could we could try the PCC now we can try the Alexander plots let me see if this works this is It's live coding, we have no time, two minutes, perfect time for live coding. I'm not sure if I have to write Germany or Deutschland, I think it's all English. And let me see if this provides, yeah, this actually works, and here is our code. But this is a really great geocoding, a free geocoding library using OpenStreetMap data to geocode our positions. So let me take the second one, and then we have, we can quickly display these two so you have an idea where it is located, and now the question is how can I get from here to here, and there is really a very simple function. I have to approximate, first thing I have to find the next node in our street network that's the nearest nodes, origin, destination, origin, destination, and now we can use this network X module and calculate shortest paths between them and let's just plot this this and we see that's the shortest path in the drive network so very easy always open sweep map I can convert that I can create GeoPanda's data frame out of it and put that also on a volume map by converting to JSON and perfectly we We have 30 seconds left. We see our street network from here to here. OK, that's basically it. This is just an announcement. If you are interested in more, there are two possibilities. I'm currently working on a multi-hour, I think it's 30 or 40 hours, online course, which is finished soon. If you are interested, you can sign up here. And I will also talk, if there is some interest, at the bar camp in Rügen. This is advertisement, so very nice opportunity to combine holidays with Python. It's actually German speaking, so that could be a problem if you don't master German. But sometimes we do have English sessions, too. OK, that's it. I'm around for questions, and thank you for your attendance.