Creating 3D Maps using Python

Speaker 1 [00:04]

Thank you for the introduction. Let me ask a question first. Who of you knows Google Maps? Everyone, of course. Who doesn't? Who used Google Maps in the past one or two weeks? Most, yeah. Who of you knows Google Earth? Who of you used Google Earth in the past one or two weeks? One, two, three. Okay. Four. Sorry. Five. Okay. Okay. I think you see the problem already. A couple of years ago, now the mouse is gone, a couple of years ago I created a virtual globe. Actually it started in 2005 and from 2011 to 2015 I developed it in WebGL. So it was a browser version. It's basically the same as Google Earth. I have an old movie here, 10 years ago, quite funny, the quality is not really the best at the moment. From today's perspective, you see some data here, can fly around, you can zoom in somewhere, you can actually jump to some predefined locations in that. It's a whole SDK, you can create such applications. It's not just a globe, it's an SDK. Okay, let's wait, it flies. it's a movie, that's not best quality. The quality was better back then. So now we go to the class here. Perfect. Really funny. A couple months ago, a student of me created a new version using another technology. That's a little bit of focus of me at the moment. I want to have this on mobile devices. Actually, it's transmitted via 5G on a mobile device and it's being rendered on a supercomputer, something like that, and you see the quality is already much better and it's on a regular smartphone. And you also see the data is much, much better than it was 10 years ago. And the cool thing actually is that all this data today is open data, I mean even 10 years or 20 years ago you had big trouble actually finding data. It sounds nice, but it isn't really. So let me talk about the goal I have at the moment. I want to create 3D models out of open GeoData. That's very important, open. We heard about OpenStreetMap before. Wonderful thing. And I really like open data. Actually, open and free is not the same. Let's not go into details there. But I want to create general 3D objects out of freely available raster. That's aerial or satellite imagery. And I want to use also a freely available elevation data for that. What I don't want to do today is I don't want to create a viewer, especially in Python. Python is a wonderful language, but I would not really create a 3D viewer in Python. There are other better ways to do that. So if you don't really know what's raster data, that's quite simple. It's basically an image. You have something like pixels and instead of pixels you do have information there. It could be just RGB or it could be elevation data. More and more cities today, they provide their data as open data and just recently, I'm from Switzerland by the way, in Switzerland the government released all base data as open data last year, actually in March 21. Also maybe thanks to the pandemic they saw how important it is to actually release data to the public because people can actually create software out of that. So we do have, in Switzerland at least, we do have these 10 cm per pixel orthophotos, so you have the whole country in this resolution. In Germany, many cities provide this data, too. Not all cities, unfortunately. For example, in Berlin, you can also get it, and in other major cities. And that's quite nice. We have 10 centimeters per pixel in every country. If you don't want to use data in Switzerland, if you want to use it in your country, maybe you're lucky. Maybe there is some data available. There are some entry points for open data. The open knowledge foundation and the data world link provides some data sets. One data set I want to look at at the moment is Landsat. It's also freely available, a satellite image data. Actually it's not only image data, there are multiple bands. We see here is just a RGB representation over the Himalayas. So what do you get from Landsat? Landsat has different versions. Now the latest is Landsat 9, but I'm still using 8 because we have many, many more images available, so we have to wait maybe one, two years, and then we can start using 9. So how it is done, you have different bands. For example, the first four bands are visible bands, so the spectrum green, blue, red, and And then you have near infrared, you have shortwave infrared, you have panchromatic, you have cirrus, and you have thermal infrared, and that's really very nice, you can do many applications with that. But now, at the moment, I'm only interested in the first three or four bands, so that's the visible spectrum, and we can actually create images out of it. For example here, that's created from Landsat, which looks quite nice, more or less, it's not really fully zoomed in, but that's freely available. However, it looks nice because it was edited with Photoshop or something like that. If you actually download it, it looks like that. So you have to actually do some work to have nice colors. It's all 16-bit imagery. And one other thing you will notice here is the black border. I made a little graphic here where you maybe see the problem if a satellite flies around this satellite is not just north aligned so to have a proper image you always get these proper aligned image always get these black borders you have to think about this as a photo done by the satellite and then of course if you align it then it is not that nice anymore. You could actually lose data and just remove these black things but if you actually do that here you will really lose too much data. You can draw some lines there and then it's all almost all gone. So that's not a good idea. So we have to download many images and, like a puzzle, put them together. There is one nice URL because Amazon provides that. Amazon provides all data on their cloud for free. You can actually freely download it. They pay for it. The first thing is if you want to use Pandas for that, actually it's nice to use Pandas for that. You can download the scene list. It's called scene list. Every image is a scene, and you can then sort by date, for example. I would actually use parse dates to acquisition date. There are several dates, processing date, acquisition date, et cetera. I parse the acquisition date, and then I take all the pictures with a cloud cover, for example, lower than 50 percent. You can also try to do three percent, one percent, zero percent, but then we will not get too many images, of course. Then you also T2 and RT, that's already pre-processed images, so with that you get all the scenes with a reduced cloud cover, and you have already a list of many, many images. The next step I did, one problem is if we take images like a puzzle, you have to imagine it's four or five years, so one image is five years old, The next is one year old. And if you put them together, it doesn't look nice because something's changed. So you should only take images which are near to each other. But that's also difficult because you want zero cloud cover. So it's a compromise you have to do. I always say it's a little bit of a manual process, too. It's not fully automated like I show it here. So at the end, I actually choose it by eye. But one thing I did is I want to have summer on all images. I don't want to have winter because the snow cover really varies among different years and then it's really not that nice, so let's use summer and for that I just take the acquisition months greater than six and lower or equal to eight. For the northern hemisphere and of course winter on the southern hemisphere is different. It doesn't really matter because, but to be perfect, I did that from December to February on the southern hemisphere. I mean, in Australia, you will not have much snow cover, so it doesn't really matter too much there, but still, I did it this way. And then we have all the scenes in the summer months. Okay. So I took two, you see. You actually don't see it because the projector is too dark. But if you see it on my screen or on a brighter screen, you will see there is actually some small difference between the two images here. And you can imagine we take more and more and more together, four, five, six, seven, to fit it together. And then after a while, you get all the coverage of the whole planet, which is really a lot of work to do that if as I said some manual looking at the images is involved there too. Maybe with machine learning that could be done something but I highly doubt it at the moment it's still the human eye looking at the images which ones fit the best. So how is it done reading raster data I recommend two libraries actually only one library if you use Python it's It's Rasterio. Rasterio is a very nice module created by Mapbox for reading geospatial rasters. So, rasters with a georeference, so every pixel has a position. And it's based on GDAL. GDAL is a geospatial data abstraction library used in every major software. If you know Qt GIS, for example, it's used there. ArcGIS, all the big GIS use GDAL. It's open source library, and it's not really Pythonic. There's a binding there, but it's really a mess from a Python perspective, let me say that. Previously, I used that in C++, that's okay, but if you have some really strange commands like open with a big O or these things or function names with big letters, that's not really common in Python, so they created more PySonic way to do that. And of course, one thing in Rust here is it's the NumPy directly NumPy, so you get NumPy arrays. That's also very nice. Actually, in GDAL you get NumPy arrays too, but it's a mess there, so don't even get started. So how is it done? It's really done in three, four lines of code. So you open dataset, you read the rasters like RGB, red, green, blue, and then you can put it together with dstack numpy and then you can display it as an image using matplotlib or whatever you like or you can store it again as another file and then you can actually access all the values there. One thing you notice maybe is the background, this black background, that's just black. It's a black pixel, zero, zero, zero. So I actually, I struggled with that a little bit because there could be a black pixel in the image. So however, it's a 16-bit image. So the probability there is really a 16-bit black pixel in an image is really low. So I create a mask out of it to just check if it's zero, 16-bit zero, and then I create this mask and then it's transparent. I will show an image soon. One thing is, of course, you have to reproject. All the images are in a different projection system. Every tile, for example, is a different UTM zone along the longitude, so you have to reproject it. I did that with GDAL tools. GDAL comes with many tools, and one tool is GDAL Warp, where you can just reproject to WGS84, for example. the EPSG code 4326. EPSG, by the way, is the European Petroleum Survey Group, and they have all the projection systems in a numeric way. So that's done for each tile. I just did it with two here, and it looks like that. I put it on a map so you can see it really fits. And the background is not black anymore, so that's nice. And the next step I didn't show here is actually to create one big image out of all the tiles and that can be done using a virtual raster data set called VRT. So let me just one part now and we're running out of time soon. Let me see, oh it's really, I feared that. So we can already retrieve images there, oh we still have 15 minutes, no problem. We can retrieve images of the whole world. If we have local data sets, as I said before, from local governments like Switzerland, Germany, they have really nice data sets. In the US, they have nice data sets too. And you can have 0.1 meter, 0.5 meter per pixel resolutions, or if you have some remote locations, then you can get the 30 meter data set from Landsat 7. And if you have a big pocket of money, you can buy data sets from satellite providers. If you say, no, I don't want to have free data sets, it's available. There are many companies providing that data, but it's really at a really high price tag. You probably saw all the images from Ukraine, for example. It's more or less, there are some providers really giving the data. They give it for free in crisis zones usually, and you could download that there, but usually you have to pay for all the data. The next step is elevation data. There are some near global data sets available. If you are looking for data sets of Antarctica, for example, very high mountains there too, that's a little bit more difficult to get. the global datasets that go up to 60, 70 degrees in latitude, so you have to specially search for them. So the two big ones are SRTM, the Shuttle Weather Topography Mission, which is already quite old, but it's a global dataset with 90 meter per pixel resolution, or the ASTRO dataset with 30 meter per pixel. And then again, you have the local datasets you can get from many countries as open data, For example, in Switzerland we have a 10 cm per pixel resolution now available which is over the whole country. Actually, not the whole country, in the mountains it's a little bit less resolution, but in the cities it's 10 cm per pixel, which is really, really nice. Actually it's wrong, it's 50 cm, but it doesn't matter. So reading elevation data, we do it the same way like we did with the RGB data sets. So I read the elevation data just with one band, and then we can access it using this numpy array and have the elevation value of every pixel. What I usually do if I have such an array is create a so-called hillshade. There is the GDAL-DEM utility, you can just call GDAL-DEM hillshade, input TIFF or whatever format you have, and you get an output TIFF with a hillshade. You can have more parameters, you can have the direction of the sun, et cetera, but it's always nice to have such a thing. It's much better than such a thing here than you see actually some of the data. So okay, now we have elevation data. Now we have image data available. Now the next step is let's create a 3D model. Creating 3D models is quite easy if you use the most easy available 3D format. In my opinion, that's the Wavefront OBJ format. It's really simple. First you have a list, starting with V, you have a list with coordinates. For example here is a cube, you have the coordinates of all edges, V the first, second, third, and then at the end you have eight vertices defined, and then after the definition of the vertices you have the faces. It's basically polygons, so you connect one, two, three, four. on this image, I started with zero. The index in the OBJ format starts with one, which is quite amusing sometimes because I forget it from time to time, and then I start with zero like this graphic here. One, two, three, four would be one, zero, one, two, three. See, it's the top polygon, and then you have the next phase, one, five, six, two, so it's zero, four, five, one, which is down there. You can have all the faces and at the end you have a 3D format. One thing you have to keep in mind is you have to be counterclockwise in the polygon definition and you can actually do that, you can create such a 3D object. Now, to create a 3D object from elevation data, you just take something like a checkerboard and then you just take the vertices and align the elevation value. I took Y as elevation, not Z. I used the X-Z axis for the plane and then Y for the elevation. Reason for that is that most local systems, or many local systems, are left-handed coordinate systems. can be adapted if you are in a country with a right-handed system. So I don't want to go into details with the code here. I want to tell you that's really simple. First, you just put the vertices there. That's the first part, this for loop. And the second for loop is the indices. So at the end, you do have 10 minutes. So let me come to an end soon. we have time for questions. So, I did that for some parts of the world. You see, this is a wireframe now just to show that it actually worked and you can do that. Then I applied some of the freely available textures, actually images, and it looks like that. Okay, so that's it, we have time for questions, but before we come to the questions, if you like Geo and you like Python, there is the yearly GeoPython conference back in reality and back in the physical world. Actually, it's also in a hybrid format, so you can enjoy it also from your sofa or wherever you want to enjoy it, so if you're interested, feel free to look more at it, GeoPython 2022. And now we have time for questions, I hope. Thank you, Martin. And we do not have any questions on Slido, so if someone here in the room has some questions, can raise their hand. For those people at home, you can use Slido. Thank you for the talk. I wanted to ask you if you use the same techniques for building, so very low scale. In that case, how do you deal with the order of vertex on triangles for internal faces? Yeah, it's basically the same. The problem is there are many city models available now, even textured. For example, Rotterdam has a wonderful big textured city model. The problem there is often the texture size is really huge, so you have to find a way to reduce the amount, and you have to work with level of detail techniques. You can't just visualize the whole city in one part. So that's quite difficult. But it's basically the same thing. You can create OBJ out of that if you have, usually it's available in ctgml and then you have to convert it for example, you can use python for that of course. And what I did a couple of years ago is we heard about OpenStreetMap, OpenStreetMap actually has some 3D information in there. So there is a roof type, for example, and that can also be really geometry. And some buildings in OpenSweepMap, they really have a 3D representation. For example, the TV tower, Alexanderplatz, is an example. It's really available fully in 3D on OpenSweepMap, and you can actually extract that and then display it on your favorite device. Yeah, I was asking just because when you have an internal, and all in a building, you have always the problem finding out which face is inside or outside, and I was wondering if there was... Yes, that's done because it's counterclockwise. You have to define all polygons counterclockwise, and then you know what's inside, what's outside. So it's back face culling, it's called the render technique. you need to have counterclockwise definition. We have a question on Slido. How long does it take to calculate all the polygons? I experience some problem with computation time for larger areas. I don't have internet here. I can just show live. It's really fast. I don't know what you did. I mean, it depends. What do you think about huge? Okay, I can't ask back, I think, but I just did it recently, a really huge error with millions of vertices, and maybe it takes a couple minutes maximum. Okay, if you say minutes is long time, okay, maybe. but for me it's not a long time to calculate. I mean, when I did the virtual globe, sometimes you had to calculate a week, two, three for a country, so that's, and it's fast. This is really parallel processing and everything, and you can be happy to be finished in a week, two. For example, if you take the whole data set of Switzerland images, that's a couple terabytes, so maybe six, seven, eight terabytes, and processing data there, processing seven terabytes, it takes time. I just did it recently in three days, and I'm happy that it only took three days. One thing, of course, if it's slow, invest in RAM. RAM is the most important thing. I do have, in my Linux machine, I have 192 gigabyte of RAM. And, no, that's, some people invest in processor power. You can buy the best processor, you can waste money on that, that's fine, but I think investing in RAM is much more important. Okay, we have other question about performance. I think you already answered. But you can still, I'm still interested in the question. I don't know if there is other people in the room, okay? Okay, and there was a question too. I don't know the answer. Oh, I don't know the answer, OK. Yeah, my question is, so for these kind of 3D models, what are the use cases? Who are the type of people that are using these things? That's a good question. As I said in the beginning, who used Google Earth in the past two weeks? And we saw it's just, OK, here's a geotrack, so it's already overrepresented in a certain way. But that's a big issue. Nobody wants 3D. And there are some use cases, for example, for solar panels and these things, okay, that's nice, but you're right, 3D is still a niche and not really ready. And I've tried to work on that already 20 years and it's still not ready. Maybe in 20 years I will give up or still try something. Does it work? Yeah. OK. Thank you for the talk. How would you compare the ease of hunting for data, for instance, to 10 years ago when you did a similar project? Is it easier to get an open or free data today? Are there databases, or do you really need to dig into the governments and whoever has what? One thing is, compared to 10 years ago, we do have much faster internet again. So data is really downloadable. I mean, 10 years ago, you had to ship hard disks, hard drives. Today, you can really download everything. You can download, I just downloaded a couple terabytes two weeks ago, so no one cared about it. If I did that 10 years ago, I would get a call from our IT department, so that changed, certainly. And also, the data availability is much higher. It's free now. For example, the data set of Switzerland costed, in the year 2005, it costed a million Euro, the data set, if you wanted to buy it. and I couldn't afford that. There was, as a researcher, as a poor researcher, we could get it for 1,000 Euro as research institute, but we couldn't provide it to other people, and that really changed. Open data has arrived. It's not perfect yet. We see some drawbacks, and we see people scared about providing data. We see conflicts in the world, and maps are still a military thing, more or less. If you know where everything is located, it has some disadvantages. Some governments don't want to provide all data sets, so that's still a fight, but more and more provided data sets, and I think that's still a good thing, and that's really changed from 10 years ago. We have finished the time. I think you can continue asking questions to Martin after the session. Thank you all. Let's give Martin another round of applause.