Introduction to spatial and spatiotemporal data in Python
Antonić, Luka
Software requirements: opengeohub/py-geo docker image (gdal, rasterio, geopandas, eumap) This tutorial explores the basics of spatial referencing, reading/writing raster datasets and vector datasets. It also shows how to access the datasets produced by the Geo-harmonizer with eumap, working with time-series datasets, and finally how to visualize the results.
 So, hello everyone, my name is Lukka, I am doing the training session about introduction to spatial and spatial temporal data in Python. So the outline of this workshop first, I would like to start with some spatial referencing basics because I found that a lot of experts that are proficient in GS actually sometimes forget the basics of coordinates, reference systems and what they mean for the data and even myself with actually background in Geodesy and it reminds myself of this sometimes. Next, we will continue with Raster IO and manipulation and after that with vector IO and manipulation in limited capacity because the rest of the workshops are mostly focused around points so we will limit ourselves to that. We will then compute a basic time series and go through some, if there is time, go through some new map convenience and performance utilities. So to continue, so what is actually spatial data? It's data in a spatial context. So when we say spatial context, we usually mean the data is placed somewhere on the surface of the earth but spatial context is more than that because we build multiple layers of abstraction around the earth surface. So we have like first there is a surface of the earth, then you have like a terrestrial reference system, then you have an ellipsoid and then after all that you have like a surface that you project ellipsoidal coordinates on. So yeah, we need all of that to put data in a spatial context. So Rasters are positioned basically with multiple parameters. You have the coordinates of the upper left angle of the Raster. You have the size of the pixels which is actually defined by four parameters, two of which are usually zero. If you see the BND, the BND, sorry. And if maybe some of you recognize this image below the left side, it's from the Wikipedia page on Esri Worldfiles which is basically what's positioning Raster in the spatial context is it's an FN transformation. And vectors are positioned in bits of a simpler way which contains more data but the data itself is written into the vectors. It's just coordinates for which vertex. So yeah, the X and Y coordinates in the context of the geospatial data are a reference to the projection surface. This is also true even for 3D data because elevation data is think of its own basically, usually. So the quality of a coordinate reference system really depends on what you are using it for. So basically for mapping, topographical mapping or web mapping, like navigation maps and that, you would use a conformal projection because it preserves angles. So things retain the shape of how you know them. But for quantifying surface properties, for example, you would use equal area projections like we use for open data science because that's what we're doing here. And why does any of this matter? Because we use data from all kinds of sources and we forget that it might not be compatible with one another. So you can have like, well, if you have continuous data, then if you project it in another CRS, then everything is probably fine because you interpolated, but if you have like, have a nominal data like, for example, land cover, you cannot interpolate it. So you do like nearest interpolation and you have pixels that are largely different shape and if you compare them with another, for example, nominal data, data set, your results would not really hold water. And also with vector data, there's things like old coordinate reference systems, like for example, in Croatia, we have before the current system, we had inherited the system that was from the Austro-Hungarian Empire, which had a local ellipsoid and people are still having trouble to this day converting data from before to our new system because you basically get vectors that are 400 meters away from where they should be and stuff like that. So yeah, all of this being said, you won't have to deal with this during the course of these training sessions because all of the data is in the same system and all of the pixels are aligned in exactly the same way. So yeah. So yeah, let's continue with the practical part. Okay so before we start, if you joined the previous session with Landra, you don't have to do this but if you are just joining us now, please update your notebooks like this. So you go to this directory and just type it pull. So okay, we can begin. So yes, again this session is very introductory. So if any of you are joining this just because you are joining the entire Python track, then sorry if you are not going to see anything useful but yeah, so let's begin. So yeah, we'll cover the basics of accessing basically and manipulating data in Python using the standard tools which are REST, RIO, Numpy, Geopandas and also UMAP which has some performance and convenience oriented functions like Landra already mentioned. So yeah, okay. Let's clear all the content. So let's begin. We will gather some files. I'm not sure if, yeah, Landra already mentioned all of this. So this is, we will use this style, this style ID, this is the data there. We can just go ahead and run this. We'll get a bunch of file paths. Then we can open the second file in the list and then we get the dataset reader. Data set reader is a REST, RIO object that enables us of course to read files and contains file metadata. So if you run the next cell, you see there's a driver, Geotiff. There's the number of bands from the file, the shape of the REST, the coordinate reference system and the transform. No dates of any as well. You can actually see all of this in the REST profile. Maybe don't print it. See the display. No, it's the same, never mind. So now if you read the data, you can put the number, the band index here. Here there's only one band, so it doesn't matter. If you do not put the band index there, it reads all the bands and just if there's only one band, it just adds one more dimension to the array. So we're going to go ahead and do this. So when you read the data, we get an umpire array. And we can do with it whatever we could with other umpire arrays. So we can do this. We can compute with the methods of the array. We can also do some stuff which we'll get to later. We will do that, compute some statistics and then plot with the plotter from your map. So yeah, this is the REST there. I think maybe the picture is a bit big. So if anyone is wondering why the fixed size is only one number, if you're familiar with Matplot Libits because the other axis is computed automatically to preserve aspect ratio of the RESTers. So yeah, here we have a, what was it? It was a blue band, I think. Yeah, blue band. Well, let's make it blue. Okay, so here's the blue band. So now we get to do some stuff with NumPy because NumPy does, enables you to do computation basically parallel under the hood. You don't have to think about parallelization and does it on its own with BLAST libraries. And one of those things we can index, we can compare arrays to a number in a vectorized manner and then index them with what we get from there. So yeah, here we'll compute the 20th and 80th percentile. If you're concerned with the usability of this, don't be, this is not about doing something particularly useful. It's just to go through some things one can do with RESTers, with RESTorio and NumPy. So yeah, we'll compute the percentiles. We'll compare the data with the high and the low percentile. Then we will combine the index with the bitwise or, so bitwise means you compare basically the indices element by element. So you have two arrays of Boolean data and you compare them element by element so you get a combined index. And then we just, yeah, we can calculate the percent of data outside the bounds by, yeah, this code. So yeah. So this is what we get. The low percentile, oh sorry, this should say 20, not 5. And this should be 80. So the low and high percentiles. And this index is, like I said, an array of Boolean values. So we can then use that array. So we will first create a copy of the data so we don't soil our original array and we will index it with this. And say that all of these pixels that are true here should be no data and then plot the results. Now again, reduce the size a little bit. Yeah, okay. So you see now that a lot of the raster is now no data. We can do something that's usually more useful and actually clip the values to the low and high percentiles like this. So we index first with the low index and say that the values should be low and then with the high index and that they should be high and we can then plot the results again. Yeah, again, let's make this blue. Yeah. Okay, actually let's not make it blue because you can see the difference here. I'm going to go back to the previous image and actually just let's plot them next to each other. Data and new data. Let's add another title here. We don't need the F there. Yeah, okay, you can see the title very well. Let's do this. Yeah, okay, so you can see the reduce of variability in the data because the right one is clipped. So there's something we can notice here as well is that this corner here and this cut out here that were previously no data are now actually clipped because no data of this raster is zero if I remember correctly, which is lower than the low percentile and we did not do anything to account for that. So what we can do is this. Restrao enables us to, we could also compare just the data to the no data value and then produce an index that way, but the more efficient way is to use the read masks from Restrao. So we read masks from the same band. We do this as type of to make it a bullion area because you get you in date, I think. If you don't do that, then we can do this again. So we get only the pixels where there's data. So we can compute only on them. So we get data only. We do again, we do the indices. We compare the data pixels to high and to low. We again clip them. We make another new data array. We fill it all of it with no data and then we can fill only the pixels that under the data mask with what is contained in this, the data only array. And again, let's actually let's plot the old data here as well again. No, sorry. No, wait. So if you can see there, no data is now accounted for. So we have made a new raster and let's now go ahead and write it into a file. So first we will create a directory with this part here. Then define the path of the output file. And now we can use Rasterio dataset writer as a context manager as one would with a regular open from the Python standard library, but aside from the output path and the mode, which is the right mode here, we need the raster profile. So if you remember the dictionary that contains all the raster metadata, well, we can just pass that here as keyboard arguments because we haven't changed basically anything other than the data itself. So the data type is the same, the CRS is the same, the raster size is the same, everything is the same except the data. So our new raster basically has the identical profile to the new one. So let's go ahead and write it. You can see we have here dataset writer object and it wrote to the file. So just to check we can use plotters with the file path to see what was written. Again we'd use the size of the plot. So there we have it. So now to some vector data. For that we will use Geopandas library that basically it's Pandas, but for spatial data it has Geopandas data frame basically differs from the Pandas data frame by having a non-optional geometry column on which you can compute geometry operations. So we will go ahead and read the land cover samples that Leandro mentioned in the session before and print some data points and then look at the points themselves. If you ignore the runtime warning there. So the CRS is the same as the raster data. Here we have the number of points and these are our points. So we will not be using any of the data in this data frame actually. In this session we will only use the geometry. So yeah it was mentioned that IPa leaflet is useful for visualizing data in spatial data in Jupyter but if you need a quick and dirty way to compare the extents of datasets this is one. So basically you can import Shapely. So box and multipolygon from Shapely geometry. Shapely is a library that Geopandas uses under the hood but used to now it's used as either that or PyGeos but yeah we won't get into that here. So you create a box from raster extents and that creates a polygon. You then do the same with the points. So you do this. Now the GascadedUnion is if you have a relatively small dataset you can make the entire data frame into a single Shapely geometry so that it visualizes easily in Jupyter and also you can compute bounds on the whole thing basically. So yeah we can do this. We have two polygons and you can see our raster is here in this little corner. And these are the points. The other points are just here and here but yeah we will just, what we will do next is take only the points that intersect with the raster. So we can do this. We create an index by using the points intersects method. So you just pass it to another geometry and you get a Boolean array that contains TrollForce depending on where the geometries intersect with what you have here. So yeah let's print the number of points we get now and let's also visualize the points subset. So yeah these are our points now. Here we can also write vector datasets. So right now we have to do this because there's bug with the current version of Fiona so if you don't do this at least we, I think with our point dataset produces an error so basically we have to do this. So what is important here other than the error wrapper is you just do this points subset to file and then you define an out path and then say which driver. You only have to define the driver if you're not writing an esri shape file because that is the default. So yeah we now wrote that to a file. So onto computing raster time series. So yeah like really Andrew already mentioned there's a helper in UMAP that enables you to read a bunch of rasters multithreadedly into a single array so you don't have to go through all the files manually and then stack the arrays and all of that. So we will go ahead and take all of the spring data set, all of the spring rasters for red and for near infrared for this style that people we were already working on. Yeah here are the shapes. As you can see these are stacked like a multi band image instead of being like this axis being on the front where you have like a list of raster datasets. So yeah we can basically just do this. Raster datasets are arrays you can compute with them directly and this is how you compute an NDI time series from the time series of red and near infrared. And we can also index one point and one pixel well and plot series for it. So yeah that's that. So we can also plot the entire series with plot rasters and I'm not sure if we have to do this anymore. I think Leandro maybe patched this last night but I'm gonna go ahead and do it anyway. Just move the axis that contains the bands to the front. So that the shape of the array looks basically like this. And yeah we are going to plot the rasters now. I think I can increase this a little bit. So yeah there we have it. We didn't plot the entire series because of the 20 rasters. We can try it but yeah. So we can now compute with the series. We can for example compute the difference of each year to the previous one. Which is basically just differencing the series from itself with the shifted index. So we take here the last 19 images and then we subtract the first 19 from them. And then we can also plot the results. So yeah now you can see differences from each year to the previous. Or you could say you want to compute differences from the first year. To this. But yeah. Let's get back to the previous year differences here. And we can use the do-map-helper function to save the rasters in parallel. So yeah we define the file names here. And we write the files. Also yeah we need to put here a path for reference raster. Which then saves rasters takes the profile from if you remember the raster profile you need when writing rasters. So instead of putting that there you just give it a reference raster to take the profile from. So yeah basically we just we unstack the data again to be a multi-band image because that's what save rasters takes. And we also should define no data here because the default is zero and the difference of Ndvi from one year to the next can definitely be zero so we should change that. We also change the data type because the GeoTIF driver does not accept float 16 which is what our data is because we computed from the start with 8-bit integer values. So it tries to keep the memory usage low by using the lowest possible float. So yeah let's save the files. And here save rasters just outputs the file names that we saved data to. So what we can do now with these files we can use the overlay helpers from map mapper to overlay all of the rasters at once with the geometries. So this optimizes things a bit and so let's run it. So is it done? It's done I think. Yeah it's done. So basically give this class a bunch of geometries as a Geopandas data frame and just a list of the files that you want intersected and it does the overlay for you. So now we have time series for each of the points. So you can see the column names are just named after the years. So let's plot the series from a single point. For a single point we are going to sort the columns here because they are not currently sourced by year. And we are going to take a series from a single point and just plot it. So yeah as you can see series at this particular point. So now onto some convenience functions from UMAP. So yeah RasterIO allows for window reading of rasters. So if you do this, if you define a window you can read raster with that window. So only the data from the window will be read. So if you do this, so you need a column and row offset and width and height in pixels. So yeah now we have a 5x5 pixel data set. But what matters here for large files is that RasterIO has to read the entire block of the raster before it reads that particular window. So if a window intersects with multiple blocks then it has to read all of the blocks before reading that window and just outputting the data. So if we look here at the number of windows that our test raster has we can just do block windows and just expand it into list because it's generator. So yeah 125 blocks. So yeah but what we can do with UMAP blocks basically we can read a larger raster by just reading it block by block window and we can not worry about that let it run in the background and just define something that we would like to do with the data from each block when it's read. So we can do this by import mapping from Shapely and take the URLs for a landsat read and near infrared for 2018. This is summer. And so we will define the file URLs and define the geometry. The geometry has to be given here in the JSON schema. So we do that with Shapely geometry mapping from our Shapely geometry. So this is the raster that we, the size of the raster that we are going to read now as opposed to 1000 by 1000 pixels that we were reading before. And this is our geometry. So we import from UMAP Parallel blocks we import the block reader, the block writer and that's basically it. All we have to do now is define NVI as a function that will be computed for each block. So we initialize the reader with the reference file and here all of the files that are read have to have the same transform and the blocks have to be identical. So this is for reading those kinds of data sets. So we give it the radio as the reference file. Then we initialize the block writer with the reader. You can also omit this and you can just initialize the writer and then that will initialize the reader with the first file. But just to illustrate because the reader is then reusable, it has to, it initializes by computing geometries for each block. So it can take a while so if you are going to reuse the reader class, the reader object rather you should initialize it on its own. So if you run that and wait a bit. So check. Oh, it's not done yet. Okay. I should run this. Yeah. Okay so if it's done. So we did, we defined our NVI function here. We defined an out file. Then we called writer write with the source paths and the destination path. And also the geometry, the function that we want to compute and should give, of course, the no data value because again I think the default is zero. The d type and now let's see if our files have written our file. Sorry. Yeah, so this is the same block but computed from very large files that are hosted in S3 buckets. So you do not have to download the entire file to work with it basically. Yeah, so the final part here, I would like to talk about some convenience abstractions around the catalog of the data produced and also used by Open Data Science. So as you know as Leander has mentioned, we have the catalog in the geo network but from a Python standpoint you would have to use CSW from OWSlib to access that directly. So here we have a catalog object that we import from EOMAP data sets. And then we can search with it. So here the catalog also has a fallback copy of the basically of the entire catalog as a CSW that's distributed with the code. So if basically if maybe you're having connection troubles and the unit per can also be a bit slow you can just use this and that's what we're going to do because it's faster and the fallback CSW is updated so there shouldn't be any difference between the two right now. So yeah, here we have a resource collection. So the output is grouped by metadata so all of the resources that have the same abstract theme and title so this is NDIA Landsat Quarterly. They're grouped basically but you can use the results object as you would list so we can just do something like this basically results. So to get the first three of them. So yeah, these are and if you look at only one of these. This is so these results objects are basically strings so you can use them as URLs for files to use with Rasterio and the UMAP helpers but they do contain additional metadata so if you want only the metadata objects you can do this. So here is the metadata. But yeah and also there is the metadata of the collection object. If you want to access the group, the metadata groups of results so if we have zero that gives us the quarterly Landsat NDVI. We index the one. We get this and oh no there's no three. Okay, this output is a bit long so we'll scroll down here. We could also use, you could exclude the results with the search so we don't have to sift through them manually later on we can just exclude with the exclude keyword. So if we exclude the trend we just basically get this. So yeah let's check results. So yeah, there's basically only the first one. Again with the long output. And we can also search by year. So here we get the quarterly Landsat for 2019. If we want to say search by multiple years so we can here input range for example from 2015 to 2020. This is non-inclusive so we should be okay. Then you have, then year is added to the metadata and you get groups of resources by year. So if we do this. So that's now just the first group. So only the rest is from 2015. So yeah we can use now these resource collections to read directly from the files. So we'll use the read rasters. And read rasters also takes a window object so we will define a window like this. So yeah the rest of the window sub module has this from bounds function which takes basically the bounds. So we will give it only the bounds that we were using so far so the extent of the test raster that's included in the VM. We will, we also have to define a transform. We can also just read from one of the files so we get the first file. We just do this. If you remember it's what the rest of the data set reader contains of all of the metadata like this. So yeah you define a window here. You just basically use read rasters as we did previously but you just say you use the spatial bin keyword argument to define a window. And this part here is basically the catalog result resource should behave like regular row strings but there's currently bug with that so we just convert all of them to strings. So we run this. And then we can plot the results. Okay this is a bit too small. Okay let's instead just plot. Well we can do that here actually. Oh sorry. Yeah so these are the windows read directly from the rasters that are hosted in S3 buckets. So okay so this is the end of the plant section so we have some time left if you have any questions if you would like to suggest that we try anything or. Are there any questions on zoom? No? Okay. Andra would you like to use the time left to go over the things from your last session that you had to skip over? I can just explain it. Okay. I think that's a good idea. Okay yeah. And yeah so you already covered some parts but I can open and explain a bit. Okay great. Thank you. So just a short recap. So we presented like how basically to access the data, read and write parts of it using the local data that we have in the virtual box like in the machine. And some of the remote data also. So this is actually the first part so to really read the points, read the raster data and in the afternoon we will use these functions but now in the context of spatial and spatial temporal machine learn. So all these functions that we present it's mainly to implement these input and output operations like deal with these GOT files. And I would like just to in my last, in my presentation the first part of this morning. So here you can see that I didn't present it but it's the same method that we use here to read these three images. You can use the same approach. You just need to change which are the file system, the file names and with these file names, these function files, it's like a utilitary function of EU map. It will sort all the files in the chronological way because the file name, you have the date inside of it. So mainly this function organize all the files in a chronological way and you can read all these files using this red raster and it's fully optimized so it reads in parallel. You can control how much threads you can use to read it and you read here all the 84 images. So here for the NIR like the infrared band of the Landsat so this is a totally harmonized product so you have it for the whole Europe and you can read it and here you have like a complete time series with four observations per each year. So we have 20 years and here it's actually 21 because it's 2000 to 2020 multiplied by four and what it's interesting so here it's just like a noon pile array with the last dimension it's the time. So each part of this array so it's like a data cube with each part of this last dimension one date but with this function here you can see I'm just plotting like a time series and I'm getting like a random pixel here and you can see the time series value so how for this specific pixel this pixel it's one position in this data cube structure you can see how it is changing. So tomorrow in the cloud optimized geative training session we will do it but with real data so we can actually access this temporal response for NDVI data directly through the cloud so it's here again you are reading the data it's just local data it's just part of Europe because we prepared it specifically for the training sessions but you can use the same approach with the same methods and you can point this to like a real data that we are hosting in the cloud and you can plot all this time series for NDVI in our app. Example and here you can compare so this is actually the structure of the data so you can compare it across different years so this method here it's just like a function to help plot the data using like the plotter the plot haster that we saw here and you can send like two years and you can compare how the image looks like so the four seasons across different years and here it's just like a operation between this this image so it's a sort of reduction so considering all this data cube you can do multiple operations and in this in this case it's just like a median you can use like maximum calculation minimum so basically it's a type of a reduction in this data structure so you have like two dimensions for the space and 84 dimensions for time each of these dimensions are specific for one season and for one year and you can derive multiple statistics about it so for example here it's you this code it actually calculates like the median value across all these images across all this time series you could have like the maximum value the minimum value so it's a easy way with the work with this data and here it's just an optimized way to access it but all the functions available in umpai you have available also to process and to deal with this data because we are creating this EU map function just to provide some helpers and to optimize for example specific parts of workflow but it's fully compatible with what we already have in python numpy scikit learn as you will see like this this afternoon so yeah this is the last part of the presentation and you can see it's possible also save the raster files and you can open it in the quantias and it's a nice way because when you save a file you need to you are reading a data put it in memory so you have a numpy array you can process and do multiple operations but to save it as image again you need to set several parameters what is the projection system and what will be the like the the data format and things like that so this method save raster it's actually used as profile one like base image and copy all these geo spatial parameters for all these images separately and you can save multiple files like that also so yeah that's that's it's mainly what we want to present this morning it's more about what data we are processing how we are reading and how we are writing and how you can access it using python so it's a preparation for the afternoon because in the afternoon we will use this data to do a proper like machine learning model and classify it and this landsat data primarily we will classify it in different land cover class so yeah I don't know if you have some questions so otherwise we can stop the session have a lunch and prepare for the afternoon sorry okay okay so let's okay let's go back here and I think if we don't if we don't have like more questions I think we could start with yeah I can actually I can execute it and because and we can yeah okay so I will open mine so I was rushing but as we have time now I will I would like to explain about this data set preparation so it's a it's a nice context for the machine learning model so mainly what we are this is mainly the data that we work until now so the landsat and R&D images so you can see here this is a nice publication about the data set but mainly we are downloading it from University of Maryland as I said and we we are getting all the like this 16 day style so for the whole world is this product provides like one integrated image every 16 days and when I say integrated image it's put together all the three landsat satellites that are available like landsat 5, 7 and 8 so it's a nice product it's a harmonized product which makes all these pixel values and when you think that there are three satellites different getting data from the earth surface of course you have different problems and it's really challenge it's really a challenge put all this data together so this is a harmonized product it's it's nice but even considering it you still have miss data and you still have like clouds so what we are doing we are getting this this data and we are organizing it for each season and mainly we have the date intervals here so for example between March 21 of March and 24 of June we will have like several images maybe six seven it depends in winter of course we will have less images because you have more problems with cloud and but it it changes across all these seasons but for this time frame for these dates you will you will have like multiple images mainly what we are doing here is we remove all the clouds so there is like a mask there is an indication of which pixels are cloud which pixels are not cloud and we are using this mask to really remove and clean the data but of course if we remove the cloud so this satellite it's not it's an optical satellite so it doesn't see across the clouds so you will have a value and a gap there so and if you have a gap if you have a miss data and no data it is a problem to multiple applications so it's it's it's better it's actually it's mandatory put some value there to use this image in this format but if you think that okay maybe for this interval I I will have like then this data will have gaps so there are two options here you can expand this window so imagine it if if I have like value pixel values here cloud here I can expand my window so I can get all the images available I don't know in 2020 for example and maybe some of these images will have like pixel value so it is not guaranteed that I will have like observations because it's a combination of when the satellite pass and if we have clouds or not what are the climate conditions but you can expand it but of course if you expand it you will have like a final image with different with pixels derived from different observations and it could mess the data so it's it's why really it's important harmonize all these observations so and but in the context of the project we decided to use we decided to use other approach so mainly we kept we kept this seasons and we get feel with you we get fill it using like different years and neighborhood seasons so and I didn't prepare the tutorial about it but we have available in the official documentation I will just show to you but we in the EU map documentation and in the official documentation we have these tutorials so and here it's like a it's exercise of how we are implementing the gap filling so this is actually the code that we run so and you can see here where is the original data and where we get filled the observations so mainly it's it use like a moving window to get observations from different seasons and from different years and fill the gaps of course you need to be aware that you are introducing like a new observation there that and in some applications it could introduce some problems so imagine it like in a case of a deforestation if you put like a gap filling maybe you are not see a real deforestation because you are using like a value from the past for example but at some point we will see that and because at some point the satellite will pass through this part of the earth to that part of the earth and you can you and you will have like a real observation for surface and you will see it so it depends of the application but here as we are mainly we are dealing with data to the past so past data since 2000 we have like some observations up forward also so what we are doing we are just using like the medium value from past observations and from future observations so for example to get filled I imagine 2020 2010 we will use observations from 2012 and from 2008 of course we will have limitations in the end of the year in the end of the time series like 2021 so and we are preparing like a publication to explain it but for the context of the project and to produce like land cover maps to the past this data it's it's it's prepared to do it and of course we have four seasons and maybe one of these seasons has actually observations without gap in the most part of the case that's that is the case actually so maybe we are gap filling for example the winter image because there are too many clouds in the season and but the summer image we have observations so it and we are using all the seasons inside of one year to produce the land cover maps so this problem will be minimized and and here you have other approach we did test three of it and one more that it's not implemented here and here you can see this is more like a linear regression between all the available observations in the cross time but you can see that our approach here it captures better the decisionality and we see these spikes because we actually we just have four observations per each year but if you think that this data it's a 30 meter data so and it was prepared for the whole Europe so we removed all the clouds and it's a nice time series and it allows to develop several applications on top of it and so let's go back and we will use a lot this data and so just to explain like the last part we did it for three percentiles so mainly inside of each season imagine you could have like four observations three observations maybe you have just one observation that's it's the case for the winter in the most part of Europe but we are calculating these three percentiles for all the seasons so it's a it's this approach it's kind of it's a straightforward way to produce to organize these different data availability in three ready to use images so for example if you have like just one pixel and for sure this three percentiles will be the same but if you have like on the Italy for example you have a better observations and a better data availability you have three different observations for each these three percentiles because you will combine all these pixels and you will produce it so the output here and of course it depends of the data availability but it's a it's a approach that allow us to generate this data and to prepare it for the machine learn and the assumption here is when we train a machine learn we will use points from the whole Europe so we will have all these different situations so we will have points with three values constant values in the three percentiles inside one season and we will have values with like multiple person with multiple values for it is different values for these three percentiles and the machine learn the model will manage it so that was the assumption here but again there are multiple ways to deal with that you can we could use like just one image for the entire year that that's a way I implemented it in the past but here we really preferred to generate the seasons and get through it and and it's it was more like a methodological decision but of course there are different strategies that could be implemented the important part and I I'm repeating it because it's not it's not these satellite images they are not like if you think in this Landsat data in the Sentinel data as they are available now they are not ready to use you need to do a lot of processing and you need to download it you need to remove the clouds you need to do some a lot of harmonization to use for example different observations from satellites so and there are there are like specific projects to do it so here we are using one of these products that already implemented it and we are building like a new kind of additional product on top of it and if you use directly the satellite images so you have a lot of work and probably during the machine learn modeling you will have you will have problems and because this even for example from the same satellite you can have like different interference different like noise from the atmosphere effects and things like that so if you are really focusing develop like models and use this data it's important to have a clean data harmonized and ready to use otherwise you need to implement that and that's the reason of why we are doing it and we are putting all this data available so to allow other people other researchers and other organizations to develop applications on top of it and of course we are improving this product so now we are generating a new version a version 1.1 and it's still this product it's not it's not perfect you can see like some artifacts still but we are really improving it and during the year so and we have the capacity to rerun it so I will as I'm talking about it I will open like one RGB image here so we can use the data portal to see it on top of my tutorial it's the easiest way to check these composites and if we go here so this is a public URL you can access and this is the RGB data and I will use the zero capacity and I will zoom in some place you can use this to yeah we can try to find one and so yeah this is nice so and we are here and you can see this is this is not like a Google or image because of course each pixel here has 30 by 30 meters so if you zoom in it it will not shift for a better image as we are as Google Earth do you have the Google that not here is the big map but it's just to show the difference but this image it's really enough to do like several environmental applications and so but what I want to show through here it's this is the summer image so it's the best mosaic that we have available because of the data availability and this is the medium value and here you are seeing actually three bands like the red green and blue and if you go back here you can see like from all the years you can see some changes and you can see still some artifacts so it's there is like rest of like some noise some artifacts because of course this processing it's not perfect we are trying to minimize the number of artifacts but sometimes there's problems in that there are problems in the cloud mask so it's really difficult remove all these problems but what it's important is we are really like creating an optimized workflow to do it and to update it and for example if we find a better way to identify the clouds we can run it and we can process again and using all these infrastructure that we are presenting here so what it's nice so we are here in Farnigin and this is Fenendal and if you see all these images you can clearly see some of the urbinary expansion here and this is only possible to see because we can have we have one image per year here and for each year we just have we have four images actually here it's just the summer but for Google Earth you have this for some areas you have multiple images but it's not really you can't do analysis with that so with this data as we are we presented today you can load the data in the memory and you can process it you can derive some statistics you can calculate like NDVI index index and for the NDVI we put it all and now I will turn it so if you look here for the NDVI you have the four seasons and I will turn on the opacity and you can do the same thing and see how the NDVI changed and as we know this now this area expanded so we had an urban expansion here and the NDVI will change because in the 2000 and the early years probably it was some kind of crop or pastry area and later it was converted to urban area so the vegetation index changed across it and if you you can check the values here so you can clearly see how it's dropping it and so probably in this specific in this specifically in this pixel maybe the conversion occurred here so and again this is a data that we prepared and we can try to identify it here yeah so you can you can see that you want considering one year to other year I will stop it here so yeah you had an expansion so we use it this data and we produce it as I presented but I will show it because it's a nice example so we produce it this land cover maps so we prepared the Landsat data we harmonized we remove at the clouds we put it available and we calculated this NDVI index and on using this data so it's not really raw data it's an harmonized product gap filled and ready to use but it's just data right it's just like we just have a spectral value like a pixel value that represents some some response of earth surface but it's not there is no like directly mean for this pixel I don't know if it's a urban area I don't know for example if it's a crop if it's a pasture and use this data we actually produce for all these years we produce it like a land cover map and we actually classify it using these spectral values we classify it as my internet here it's not helping you see yeah now I have it so here it's like this is the result of the machine learn model it's so you can see that this area were like it was a pastry I think and so what it's interesting here if you look the image you can identify it but yeah this because we all most like the most part knows like what is a urban area we can see that maybe the we can see the building growing up across the landscape but for a machine learn model and to produce like a product on top of it we really need to train a model we need to provide training samples so there is a whole workflow on top of this data and if you put here so here you can have like you can see all the land cover class and when it was converted to the urban area so and this this is like the result of the land cover classification we will do it on the afternoon so we will produce something like this as output for a specific for the tiles that we we have and using the same the same workflow and you have other classes here and yeah that's it's a more like it's a advertisement a way to try to stimulate you for the next sessions until the afternoon and we will really process this data and produce one like the rival product on top of this lands up data and you can you can use it for other applications also so inside of the project we are using this data to map where are the three species for example you can use this data to identify for example and to map the soil properties and you can there are several applications to build on it but it's important to emphasize that to depending of the scale depending of the area of your interest if you plan to do something like that for a whole country for example all this data preparation all this data immunization it's it's hard working and it's it's hard work and we are providing it as service for different users and different organizations and possible applications okay we have more questions or any questions actually okay okay it's about the difference between the end DVI why it's subtracting the previous year here so it's the first year yes I think it's related with yours I can explain yep okay it's more related with the Python syntax okay so so what I'm reading here is it looks to you like we're subtracting the last year from the first year right because yeah this the Python syntax is because these this means that so we're not using the first year we're using every year after the first because the first year has the index zero and the second one has the index one and then if you say one and then then this the yeah yeah okay you get it now great so was that it oh yeah yeah indexing from one is is better but yeah I'm spiteful in Harris from a long line of of working on limited hardware and all of that so yeah no problem yeah yeah indexes start from one it's weird for me but it's how our manager think okay so more questions are let's stop and now I think it's time to the lunch okay thank you all and see afternoon we see each other in afternoon
 So, hello everyone, my name is Lukka, I am doing the training session about introduction  to spatial and spatial temporal data in Python.  So the outline of this workshop first, I would like to start with some spatial referencing  basics because I found that a lot of experts that are proficient in GS actually sometimes  forget the basics of coordinates, reference systems and what they mean for the data and  even myself with actually background in Geodesy and it reminds myself of this sometimes.  Next, we will continue with Raster IO and manipulation and after that with vector IO and manipulation  in limited capacity because the rest of the workshops are mostly focused around points  so we will limit ourselves to that.  We will then compute a basic time series and go through some, if there is time, go through  some new map convenience and performance utilities.  So to continue, so what is actually spatial data?  It's data in a spatial context.  So when we say spatial context, we usually mean the data is placed somewhere on the surface  of the earth but spatial context is more than that because we build multiple layers of abstraction  around the earth surface.  So we have like first there is a surface of the earth, then you have like a terrestrial  reference system, then you have an ellipsoid and then after all that you have like a surface  that you project ellipsoidal coordinates on.  So yeah, we need all of that to put data in a spatial context.  So Rasters are positioned basically with multiple parameters.  You have the coordinates of the upper left angle of the Raster.  You have the size of the pixels which is actually defined by four parameters, two of which are  usually zero.  If you see the BND, the BND, sorry.  And if maybe some of you recognize this image below the left side, it's from the Wikipedia  page on Esri Worldfiles which is basically what's positioning Raster in the spatial context  is it's an FN transformation.  And vectors are positioned in bits of a simpler way which contains more data but the data  itself is written into the vectors.  It's just coordinates for which vertex.  So yeah, the X and Y coordinates in the context of the geospatial data are a reference to  the projection surface.  This is also true even for 3D data because elevation data is think of its own basically,  usually.  So the quality of a coordinate reference system really depends on what you are using it for.  So basically for mapping, topographical mapping or web mapping, like navigation maps and that,  you would use a conformal projection because it preserves angles.  So things retain the shape of how you know them.  But for quantifying surface properties, for example, you would use equal area projections  like we use for open data science because that's what we're doing here.  And why does any of this matter?  Because we use data from all kinds of sources and we forget that it might not be compatible  with one another.  So you can have like, well, if you have continuous data, then if you project it in another CRS,  then everything is probably fine because you interpolated, but if you have like, have a  nominal data like, for example, land cover, you cannot interpolate it.  So you do like nearest interpolation and you have pixels that are largely different shape  and if you compare them with another, for example, nominal data, data set, your results  would not really hold water.  And also with vector data, there's things like old coordinate reference systems, like  for example, in Croatia, we have before the current system, we had inherited the system  that was from the Austro-Hungarian Empire, which had a local ellipsoid and people are  still having trouble to this day converting data from before to our new system because  you basically get vectors that are 400 meters away from where they should be and stuff like  that.  So yeah, all of this being said, you won't have to deal with this during the course of  these training sessions because all of the data is in the same system and all of the pixels  are aligned in exactly the same way.  So yeah.  So yeah, let's continue with the practical part.  Okay so before we start, if you joined the previous session with Landra, you don't have  to do this but if you are just joining us now, please update your notebooks like this.  So you go to this directory and just type it pull.  So okay, we can begin.  So yes, again this session is very introductory.  So if any of you are joining this just because you are joining the entire Python track, then  sorry if you are not going to see anything useful but yeah, so let's begin.  So yeah, we'll cover the basics of accessing basically and manipulating data in Python using  the standard tools which are REST, RIO, Numpy, Geopandas and also UMAP which has some performance  and convenience oriented functions like Landra already mentioned.  So yeah, okay.  Let's clear all the content.  So let's begin.  We will gather some files.  I'm not sure if, yeah, Landra already mentioned all of this.  So this is, we will use this style, this style ID, this is the data there.  We can just go ahead and run this.  We'll get a bunch of file paths.  Then we can open the second file in the list and then we get the dataset reader.  Data set reader is a REST, RIO object that enables us of course to read files and contains  file metadata.  So if you run the next cell, you see there's a driver, Geotiff.  There's the number of bands from the file, the shape of the REST, the coordinate reference  system and the transform.  No dates of any as well.  You can actually see all of this in the REST profile.  Maybe don't print it.  See the display.  No, it's the same, never mind.  So now if you read the data, you can put the number, the band index here.  Here there's only one band, so it doesn't matter.  If you do not put the band index there, it reads all the bands and just if there's only  one band, it just adds one more dimension to the array.  So we're going to go ahead and do this.  So when you read the data, we get an umpire array.  And we can do with it whatever we could with other umpire arrays.  So we can do this.  We can compute with the methods of the array.  We can also do some stuff which we'll get to later.  We will do that, compute some statistics and then plot with the plotter from your map.  So yeah, this is the REST there.  I think maybe the picture is a bit big.  So if anyone is wondering why the fixed size is only one number, if you're familiar with  Matplot Libits because the other axis is computed automatically to preserve aspect ratio of  the RESTers.  So yeah, here we have a, what was it?  It was a blue band, I think.  Yeah, blue band.  Well, let's make it blue.  Okay, so here's the blue band.  So now we get to do some stuff with NumPy because NumPy does, enables you to do computation  basically parallel under the hood.  You don't have to think about parallelization and does it on its own with BLAST libraries.  And one of those things we can index, we can compare arrays to a number in a vectorized  manner and then index them with what we get from there.  So yeah, here we'll compute the 20th and 80th percentile.  If you're concerned with the usability of this, don't be, this is not about doing something  particularly useful.  It's just to go through some things one can do with RESTers, with RESTorio and NumPy.  So yeah, we'll compute the percentiles.  We'll compare the data with the high and the low percentile.  Then we will combine the index with the bitwise or, so bitwise means you compare basically  the indices element by element.  So you have two arrays of Boolean data and you compare them element by element so you  get a combined index.  And then we just, yeah, we can calculate the percent of data outside the bounds by, yeah,  this code.  So yeah.  So this is what we get.  The low percentile, oh sorry, this should say 20, not 5.  And this should be 80.  So the low and high percentiles.  And this index is, like I said, an array of Boolean values.  So we can then use that array.  So we will first create a copy of the data so we don't soil our original array and we  will index it with this.  And say that all of these pixels that are true here should be no data and then plot the  results.  Now again, reduce the size a little bit.  Yeah, okay.  So you see now that a lot of the raster is now no data.  We can do something that's usually more useful and actually clip the values to the low and  high percentiles like this.  So we index first with the low index and say that the values should be low and then with  the high index and that they should be high and we can then plot the results again.  Yeah, again, let's make this blue.  Yeah.  Okay, actually let's not make it blue because you can see the difference here.  I'm going to go back to the previous image and actually just let's plot them next to  each other.  Data and new data.  Let's add another title here.  We don't need the F there.  Yeah, okay, you can see the title very well.  Let's do this.  Yeah, okay, so you can see the reduce of variability in the data because the right  one is clipped.  So there's something we can notice here as well is that this corner here and this cut  out here that were previously no data are now actually clipped because no data of this  raster is zero if I remember correctly, which is lower than the low percentile and we did  not do anything to account for that.  So what we can do is this.  Restrao enables us to, we could also compare just the data to the no data value and then  produce an index that way, but the more efficient way is to use the read masks from Restrao.  So we read masks from the same band.  We do this as type of to make it a bullion area because you get you in date, I think.  If you don't do that, then we can do this again.  So we get only the pixels where there's data.  So we can compute only on them.  So we get data only.  We do again, we do the indices.  We compare the data pixels to high and to low.  We again clip them.  We make another new data array.  We fill it all of it with no data and then we can fill only the pixels that under the  data mask with what is contained in this, the data only array.  And again, let's actually let's plot the old data here as well again.  No, sorry.  No, wait.  So if you can see there, no data is now accounted for.  So we have made a new raster and let's now go ahead and write it into a file.  So first we will create a directory with this part here.  Then define the path of the output file.  And now we can use Rasterio dataset writer as a context manager as one would with a regular  open from the Python standard library, but aside from the output path and the mode, which  is the right mode here, we need the raster profile.  So if you remember the dictionary that contains all the raster metadata, well, we can just  pass that here as keyboard arguments because we haven't changed basically anything other  than the data itself.  So the data type is the same, the CRS is the same, the raster size is the same, everything  is the same except the data.  So our new raster basically has the identical profile to the new one.  So let's go ahead and write it.  You can see we have here dataset writer object and it wrote to the file.  So just to check we can use plotters with the file path to see what was written.  Again we'd use the size of the plot.  So there we have it.  So now to some vector data.  For that we will use Geopandas library that basically it's Pandas, but for spatial data  it has Geopandas data frame basically differs from the Pandas data frame by having a non-optional  geometry column on which you can compute geometry operations.  So we will go ahead and read the land cover samples that Leandro mentioned in the session  before and print some data points and then look at the points themselves.  If you ignore the runtime warning there.  So the CRS is the same as the raster data.  Here we have the number of points and these are our points.  So we will not be using any of the data in this data frame actually.  In this session we will only use the geometry.  So yeah it was mentioned that IPa leaflet is useful for visualizing data in spatial  data in Jupyter but if you need a quick and dirty way to compare the extents of datasets  this is one.  So basically you can import Shapely.  So box and multipolygon from Shapely geometry.  Shapely is a library that Geopandas uses under the hood but used to now it's used as either  that or PyGeos but yeah we won't get into that here.  So you create a box from raster extents and that creates a polygon.  You then do the same with the points.  So you do this.  Now the GascadedUnion is if you have a relatively small dataset you can make the entire data  frame into a single Shapely geometry so that it visualizes easily in Jupyter and also you  can compute bounds on the whole thing basically.  So yeah we can do this.  We have two polygons and you can see our raster is here in this little corner.  And these are the points.  The other points are just here and here but yeah we will just, what we will do next is  take only the points that intersect with the raster.  So we can do this.  We create an index by using the points intersects method.  So you just pass it to another geometry and you get a Boolean array that contains TrollForce  depending on where the geometries intersect with what you have here.  So yeah let's print the number of points we get now and let's also visualize the points  subset.  So yeah these are our points now.  Here we can also write vector datasets.  So right now we have to do this because there's bug with the current version of Fiona so if  you don't do this at least we, I think with our point dataset produces an error so basically  we have to do this.  So what is important here other than the error wrapper is you just do this points subset to  file and then you define an out path and then say which driver.  You only have to define the driver if you're not writing an esri shape file because that  is the default.  So yeah we now wrote that to a file.  So onto computing raster time series.  So yeah like really Andrew already mentioned there's a helper in UMAP that enables you  to read a bunch of rasters multithreadedly into a single array so you don't have to go  through all the files manually and then stack the arrays and all of that.  So we will go ahead and take all of the spring data set, all of the spring rasters for red  and for near infrared for this style that people we were already working on.  Yeah here are the shapes.  As you can see these are stacked like a multi band image instead of being like this axis  being on the front where you have like a list of raster datasets.  So yeah we can basically just do this.  Raster datasets are arrays you can compute with them directly and this is how you compute  an NDI time series from the time series of red and near infrared.  And we can also index one point and one pixel well and plot series for it.  So yeah that's that.  So we can also plot the entire series with plot rasters and I'm not sure if we have to  do this anymore.  I think Leandro maybe patched this last night but I'm gonna go ahead and do it anyway.  Just move the axis that contains the bands to the front.  So that the shape of the array looks basically like this.  And yeah we are going to plot the rasters now.  I think I can increase this a little bit.  So yeah there we have it.  We didn't plot the entire series because of the 20 rasters.  We can try it but yeah.  So we can now compute with the series.  We can for example compute the difference of each year to the previous one.  Which is basically just differencing the series from itself with the shifted index.  So we take here the last 19 images and then we subtract the first 19 from them.  And then we can also plot the results.  So yeah now you can see differences from each year to the previous.  Or you could say you want to compute differences from the first year.  To this.  But yeah.  Let's get back to the previous year differences here.  And we can use the do-map-helper function to save the rasters in parallel.  So yeah we define the file names here.  And we write the files.  Also yeah we need to put here a path for reference raster.  Which then saves rasters takes the profile from if you remember the raster profile you  need when writing rasters.  So instead of putting that there you just give it a reference raster to take the profile  from.  So yeah basically we just we unstack the data again to be a multi-band image because that's  what save rasters takes.  And we also should define no data here because the default is zero and the difference of  Ndvi from one year to the next can definitely be zero so we should change that.  We also change the data type because the GeoTIF driver does not accept float 16 which is what  our data is because we computed from the start with 8-bit integer values.  So it tries to keep the memory usage low by using the lowest possible float.  So yeah let's save the files.  And here save rasters just outputs the file names that we saved data to.  So what we can do now with these files we can use the overlay helpers from map mapper  to overlay all of the rasters at once with the geometries.  So this optimizes things a bit and so let's run it.  So is it done?  It's done I think.  Yeah it's done.  So basically give this class a bunch of geometries as a Geopandas data frame and just a list  of the files that you want intersected and it does the overlay for you.  So now we have time series for each of the points.  So you can see the column names are just named after the years.  So let's plot the series from a single point.  For a single point we are going to sort the columns here because they are not currently  sourced by year.  And we are going to take a series from a single point and just plot it.  So yeah as you can see series at this particular point.  So now onto some convenience functions from UMAP.  So yeah RasterIO allows for window reading of rasters.  So if you do this, if you define a window you can read raster with that window.  So only the data from the window will be read.  So if you do this, so you need a column and row offset and width and height in pixels.  So yeah now we have a 5x5 pixel data set.  But what matters here for large files is that RasterIO has to read the entire block of the  raster before it reads that particular window.  So if a window intersects with multiple blocks then it has to read all of the blocks before  reading that window and just outputting the data.  So if we look here at the number of windows that our test raster has we can just do block  windows and just expand it into list because it's generator.  So yeah 125 blocks.  So yeah but what we can do with UMAP blocks basically we can read a larger raster by just  reading it block by block window and we can not worry about that let it run in the background  and just define something that we would like to do with the data from each block when it's  read.  So we can do this by import mapping from Shapely and take the URLs for a landsat read and near  infrared for 2018.  This is summer.  And so we will define the file URLs and define the geometry.  The geometry has to be given here in the JSON schema.  So we do that with Shapely geometry mapping from our Shapely geometry.  So this is the raster that we, the size of the raster that we are going to read now as  opposed to 1000 by 1000 pixels that we were reading before.  And this is our geometry.  So we import from UMAP Parallel blocks we import the block reader, the block writer and  that's basically it.  All we have to do now is define NVI as a function that will be computed for each block.  So we initialize the reader with the reference file and here all of the files that are read  have to have the same transform and the blocks have to be identical.  So this is for reading those kinds of data sets.  So we give it the radio as the reference file.  Then we initialize the block writer with the reader.  You can also omit this and you can just initialize the writer and then that will initialize the  reader with the first file.  But just to illustrate because the reader is then reusable, it has to, it initializes  by computing geometries for each block.  So it can take a while so if you are going to reuse the reader class, the reader object  rather you should initialize it on its own.  So if you run that and wait a bit.  So check.  Oh, it's not done yet.  Okay.  I should run this.  Yeah.  Okay so if it's done.  So we did, we defined our NVI function here.  We defined an out file.  Then we called writer write with the source paths and the destination path.  And also the geometry, the function that we want to compute and should give, of course,  the no data value because again I think the default is zero.  The d type and now let's see if our files have written our file.  Sorry.  Yeah, so this is the same block but computed from very large files that are hosted in S3  buckets.  So you do not have to download the entire file to work with it basically.  Yeah, so the final part here, I would like to talk about some convenience abstractions  around the catalog of the data produced and also used by Open Data Science.  So as you know as Leander has mentioned, we have the catalog in the geo network but from  a Python standpoint you would have to use CSW from OWSlib to access that directly.  So here we have a catalog object that we import from EOMAP data sets.  And then we can search with it.  So here the catalog also has a fallback copy of the basically of the entire catalog as  a CSW that's distributed with the code.  So if basically if maybe you're having connection troubles and the unit per can also be a bit  slow you can just use this and that's what we're going to do because it's faster and  the fallback CSW is updated so there shouldn't be any difference between the two right now.  So yeah, here we have a resource collection.  So the output is grouped by metadata so all of the resources that have the same abstract  theme and title so this is NDIA Landsat Quarterly.  They're grouped basically but you can use the results object as you would list so we  can just do something like this basically results.  So to get the first three of them.  So yeah, these are and if you look at only one of these.  This is so these results objects are basically strings so you can use them as URLs for files  to use with Rasterio and the UMAP helpers but they do contain additional metadata so  if you want only the metadata objects you can do this.  So here is the metadata.  But yeah and also there is the metadata of the collection object.  If you want to access the group, the metadata groups of results so if we have zero that  gives us the quarterly Landsat NDVI.  We index the one.  We get this and oh no there's no three.  Okay, this output is a bit long so we'll scroll down here.  We could also use, you could exclude the results with the search so we don't have to  sift through them manually later on we can just exclude with the exclude keyword.  So if we exclude the trend we just basically get this.  So yeah let's check results.  So yeah, there's basically only the first one.  Again with the long output.  And we can also search by year.  So here we get the quarterly Landsat for 2019.  If we want to say search by multiple years so we can here input range for example from  2015 to 2020.  This is non-inclusive so we should be okay.  Then you have, then year is added to the metadata and you get groups of resources by year.  So if we do this.  So that's now just the first group.  So only the rest is from 2015.  So yeah we can use now these resource collections to read directly from the files.  So we'll use the read rasters.  And read rasters also takes a window object so we will define a window like this.  So yeah the rest of the window sub module has this from bounds function which takes  basically the bounds.  So we will give it only the bounds that we were using so far so the extent of the test  raster that's included in the VM.  We will, we also have to define a transform.  We can also just read from one of the files so we get the first file.  We just do this.  If you remember it's what the rest of the data set reader contains of all of the metadata  like this.  So yeah you define a window here.  You just basically use read rasters as we did previously but you just say you use the  spatial bin keyword argument to define a window.  And this part here is basically the catalog result resource should behave like regular  row strings but there's currently bug with that so we just convert all of them to strings.  So we run this.  And then we can plot the results.  Okay this is a bit too small.  Okay let's instead just plot.  Well we can do that here actually.  Oh sorry.  Yeah so these are the windows read directly from the rasters that are hosted in S3 buckets.  So okay so this is the end of the plant section so we have some time left if you have any  questions if you would like to suggest that we try anything or.  Are there any questions on zoom?  No?  Okay.  Andra would you like to use the time left to go over the things from your last session  that you had to skip over?  I can just explain it.  Okay.  I think that's a good idea.  Okay yeah.  And yeah so you already covered some parts but I can open and explain a bit.  Okay great.  Thank you.  So just a short recap.  So we presented like how basically to access the data, read and write parts of it using  the local data that we have in the virtual box like in the machine.  And some of the remote data also.  So this is actually the first part so to really read the points, read the raster data and  in the afternoon we will use these functions but now in the context of spatial and spatial  temporal machine learn.  So all these functions that we present it's mainly to implement these input and output  operations like deal with these GOT files.  And I would like just to in my last, in my presentation the first part of this morning.  So here you can see that I didn't present it but it's the same method that we use here  to read these three images.  You can use the same approach.  You just need to change which are the file system, the file names and with these file  names, these function files, it's like a utilitary function of EU map.  It will sort all the files in the chronological way because the file name, you have the date  inside of it.  So mainly this function organize all the files in a chronological way and you can read all  these files using this red raster and it's fully optimized so it reads in parallel.  You can control how much threads you can use to read it and you read here all the 84 images.  So here for the NIR like the infrared band of the Landsat so this is a totally harmonized  product so you have it for the whole Europe and you can read it and here you have like  a complete time series with four observations per each year.  So we have 20 years and here it's actually 21 because it's 2000 to 2020 multiplied by  four and what it's interesting so here it's just like a noon pile array with the last  dimension it's the time.  So each part of this array so it's like a data cube with each part of this last dimension  one date but with this function here you can see I'm just plotting like a time series and  I'm getting like a random pixel here and you can see the time series value so how for this  specific pixel this pixel it's one position in this data cube structure you can see how  it is changing.  So tomorrow in the cloud optimized geative training session we will do it but with real  data so we can actually access this temporal response for NDVI data directly through the  cloud so it's here again you are reading the data it's just local data it's just part of  Europe because we prepared it specifically for the training sessions but you can use  the same approach with the same methods and you can point this to like a real data that  we are hosting in the cloud and you can plot all this time series for NDVI in our app.  Example and here you can compare so this is actually the structure of the data so you  can compare it across different years so this method here it's just like a function to help  plot the data using like the plotter the plot haster that we saw here and you can send  like two years and you can compare how the image looks like so the four seasons across  different years and here it's just like a operation between this this image so it's  a sort of reduction so considering all this data cube you can do multiple operations and  in this in this case it's just like a median you can use like maximum calculation minimum  so basically it's a type of a reduction in this data structure so you have like two dimensions  for the space and 84 dimensions for time each of these dimensions are specific for one season  and for one year and you can derive multiple statistics about it so for example here it's  you this code it actually calculates like the median value across all these images across  all this time series you could have like the maximum value the minimum value so it's a  easy way with the work with this data and here it's just an optimized way to access  it but all the functions available in umpai you have available also to process and to  deal with this data because we are creating this EU map function just to provide some  helpers and to optimize for example specific parts of workflow but it's fully compatible  with what we already have in python numpy scikit learn as you will see like this this afternoon  so yeah this is the last part of the presentation and you can see it's possible also save the  raster files and you can open it in the quantias and it's a nice way because when you save  a file you need to you are reading a data put it in memory so you have a numpy array  you can process and do multiple operations but to save it as image again you need to  set several parameters what is the projection system and what will be the like the the data  format and things like that so this method save raster it's actually used as profile  one like base image and copy all these geo spatial parameters for all these images separately  and you can save multiple files like that also so yeah that's that's it's mainly what  we want to present this morning it's more about what data we are processing how we are  reading and how we are writing and how you can access it using python so it's a preparation  for the afternoon because in the afternoon we will use this data to do a proper like  machine learning model and classify it and this landsat data primarily we will classify  it in different land cover class so yeah I don't know if you have some questions so  otherwise we can stop the session have a lunch and prepare for the afternoon sorry okay okay  so let's okay let's go back here and I think if we don't if we don't have like more questions  I think we could start with yeah I can actually I can execute it and because and we can yeah  okay so I will open mine so I was rushing but as we have time now I will I would like  to explain about this data set preparation so it's a it's a nice context for the machine  learning model so mainly what we are this is mainly the data that we work until now so  the landsat and R&D images so you can see here this is a nice publication about the  data set but mainly we are downloading it from University of Maryland as I said and  we we are getting all the like this 16 day style so for the whole world is this product  provides like one integrated image every 16 days and when I say integrated image it's  put together all the three landsat satellites that are available like landsat 5, 7 and 8  so it's a nice product it's a harmonized product which makes all these pixel values  and when you think that there are three satellites different getting data from the earth surface  of course you have different problems and it's really challenge it's really a challenge  put all this data together so this is a harmonized product it's it's nice but even considering  it you still have miss data and you still have like clouds so what we are doing we are  getting this this data and we are organizing it for each season and mainly we have the  date intervals here so for example between March 21 of March and 24 of June we will have  like several images maybe six seven it depends in winter of course we will have less images  because you have more problems with cloud and but it it changes across all these seasons  but for this time frame for these dates you will you will have like multiple images mainly  what we are doing here is we remove all the clouds so there is like a mask there is an  indication of which pixels are cloud which pixels are not cloud and we are using this  mask to really remove and clean the data but of course if we remove the cloud so this satellite  it's not it's an optical satellite so it doesn't see across the clouds so you will have a value  and a gap there so and if you have a gap if you have a miss data and no data it is a problem  to multiple applications so it's it's it's better it's actually it's mandatory put some  value there to use this image in this format but if you think that okay maybe for this  interval I I will have like then this data will have gaps so there are two options here  you can expand this window so imagine it if if I have like value pixel values here cloud  here I can expand my window so I can get all the images available I don't know in 2020  for example and maybe some of these images will have like pixel value so it is not guaranteed  that I will have like observations because it's a combination of when the satellite  pass and if we have clouds or not what are the climate conditions but you can expand  it but of course if you expand it you will have like a final image with different with  pixels derived from different observations and it could mess the data so it's it's why  really it's important harmonize all these observations so and but in the context of the project  we decided to use we decided to use other approach so mainly we kept we kept this seasons  and we get feel with you we get fill it using like different years and neighborhood seasons so and I  didn't prepare the tutorial about it but we have available in the official documentation  I will just show to you but we in the EU map documentation and in the official documentation  we have these tutorials so and here it's like a it's exercise of how we are implementing  the gap filling so this is actually the code that we run so and you can see here where is  the original data and where we get filled the observations so mainly it's it use like a moving  window to get observations from different seasons and from different years and fill the gaps of  course you need to be aware that you are introducing like a new observation there that and in some  applications it could introduce some problems so imagine it like in a case of a deforestation if  you put like a gap filling maybe you are not see a real deforestation because you are using like a  value from the past for example but at some point we will see that and because at some point the  satellite will pass through this part of the earth to that part of the earth and you can you and you  will have like a real observation for surface and you will see it so it depends of the application  but here as we are mainly we are dealing with data to the past so past data since 2000 we have  like some observations up forward also so what we are doing we are just using like the medium  value from past observations and from future observations so for example to get filled I  imagine 2020 2010 we will use observations from 2012 and from 2008 of course we will have limitations  in the end of the year in the end of the time series like 2021 so and we are preparing like  a publication to explain it but for the context of the project and to produce like land cover maps  to the past this data it's it's it's prepared to do it and of course we have four seasons and maybe  one of these seasons has actually observations without gap in the most part of the case that's  that is the case actually so maybe we are gap filling for example the winter image because  there are too many clouds in the season and but the summer image we have observations so it and we  are using all the seasons inside of one year to produce the land cover maps so this problem will  be minimized and and here you have other approach we did test three of it and one more that it's  not implemented here and here you can see this is more like a linear regression between all the  available observations in the cross time but you can see that our approach here it captures better  the decisionality and we see these spikes because we actually we just have four observations per  each year but if you think that this data it's a 30 meter data so and it was prepared for the whole  Europe so we removed all the clouds and it's a nice time series and it allows to develop several  applications on top of it and so let's go back and we will use a lot this data and so just to  explain like the last part we did it for three percentiles so mainly inside of each season  imagine you could have like four observations three observations maybe you have just one observation  that's it's the case for the winter in the most part of Europe but we are calculating these three  percentiles for all the seasons so it's a it's this approach it's kind of it's a straightforward way  to produce to organize these different data availability in three ready to use images so for  example if you have like just one pixel and for sure this three percentiles will be the same but  if you have like on the Italy for example you have a better observations and a better data  availability you have three different observations for each these three percentiles because you will  combine all these pixels and you will produce it so the output here and of course it depends  of the data availability but it's a it's a approach that allow us to generate this data and to prepare  it for the machine learn and the assumption here is when we train a machine learn we will use points  from the whole Europe so we will have all these different situations so we will have points with  three values constant values in the three percentiles inside one season and we will have values with  like multiple person with multiple values for it is different values for these three percentiles  and the machine learn the model will manage it so that was the assumption here but again there are  multiple ways to deal with that you can we could use like just one image for the entire year that  that's a way I implemented it in the past but here we really preferred to generate the seasons and  get through it and and it's it was more like a methodological decision but of course there are  different strategies that could be implemented the important part and I I'm repeating it because  it's not it's not these satellite images they are not like if you think in this Landsat data in  the Sentinel data as they are available now they are not ready to use you need to do a lot of  processing and you need to download it you need to remove the clouds you need to do some a lot of  harmonization to use for example different observations from satellites so and there are  there are like specific projects to do it so here we are using one of these products that already  implemented it and we are building like a new kind of additional product on top of it and if you  use directly the satellite images so you have a lot of work and probably during the machine  learn modeling you will have you will have problems and because this even for example from the same  satellite you can have like different interference different like noise from the atmosphere effects  and things like that so if you are really focusing develop like models and use this data it's important  to have a clean data harmonized and ready to use otherwise you need to implement that and that's  the reason of why we are doing it and we are putting all this data available so to allow  other people other researchers and other organizations to develop applications on top  of it and of course we are improving this product so now we are generating a new version a version  1.1 and it's still this product it's not it's not perfect you can see like some artifacts still but  we are really improving it and during the year so and we have the capacity to rerun it so I will  as I'm talking about it I will open like one RGB image here so we can use the data portal to see it  on top of my tutorial it's the easiest way to check these composites and if we go here so this is a  public URL you can access and this is the RGB data and I will use the zero capacity and I will zoom  in some place you can use this to yeah we can try to find one and so yeah this is nice so and we are  here and you can see this is this is not like a Google or image because of course each pixel here has  30 by 30 meters so if you zoom in it it will not shift for a better image as we are as Google Earth  do you have the Google that not here is the big map but it's just to show the difference but this  image it's really enough to do like several environmental applications and so but what I  want to show through here it's this is the summer image so it's the best mosaic that we have available  because of the data availability and this is the medium value and here you are seeing actually three  bands like the red green and blue and if you go back here you can see like from all the years you  can see some changes and you can see still some artifacts so it's there is like rest of like some  noise some artifacts because of course this processing it's not perfect we are trying to  minimize the number of artifacts but sometimes there's problems in that there are problems in  the cloud mask so it's really difficult remove all these problems but what it's important is we are  really like creating an optimized workflow to do it and to update it and for example if we find a  better way to identify the clouds we can run it and we can process again and using all these  infrastructure that we are presenting here so what it's nice so we are here in Farnigin and this  is Fenendal and if you see all these images you can clearly see some of the urbinary expansion  here and this is only possible to see because we can have we have one image per year here and for  each year we just have we have four images actually here it's just the summer but for Google Earth  you have this for some areas you have multiple images but it's not really you can't do analysis  with that so with this data as we are we presented today you can load the data in the memory and you  can process it you can derive some statistics you can calculate like NDVI index index and  for the NDVI we put it all and now I will turn it so if you look here for the NDVI you have the  four seasons and I will turn on the opacity and you can do the same thing and see how the NDVI  changed and as we know this now this area expanded so we had an urban expansion here and the NDVI  will change because in the 2000 and the early years probably it was some kind of crop or  pastry area and later it was converted to urban area so the vegetation index changed  across it and if you you can check the values here so you can clearly see how it's dropping it and  so probably in this specific in this specifically in this pixel maybe the conversion occurred here  so and again this is a data that we prepared and we can try to identify it here yeah so you can  you can see that you want considering one year to other year I will stop it here so yeah  you had an expansion so we use it this data and we produce it as I presented but I will show it  because it's a nice example so we produce it this land cover maps so we prepared the Landsat  data we harmonized we remove at the clouds we put it available and we calculated this NDVI index  and on using this data so it's not really raw data it's an harmonized product gap filled and ready  to use but it's just data right it's just like we just have a spectral value like a pixel value  that represents some some response of earth surface but it's not there is no like directly mean for  this pixel I don't know if it's a urban area I don't know for example if it's a crop if it's a  pasture and use this data we actually produce for all these years we produce it like a land cover  map and we actually classify it using these spectral values we classify it as my internet here it's  not helping you see yeah now I have it so here it's like this is the result of the machine  learn model it's so you can see that this area were like it was a pastry I think and so what it's  interesting here if you look the image you can identify it but yeah this because we all most  like the most part knows like what is a urban area we can see that maybe the we can see the  building growing up across the landscape but for a machine learn model and to produce like a product  on top of it we really need to train a model we need to provide training samples so there is a  whole workflow on top of this data and if you put here so here you can have like you can see all  the land cover class and when it was converted to the urban area so and this this is like the  result of the land cover classification we will do it on the afternoon so we will produce something  like this as output for a specific for the tiles that we we have and using the same the same workflow  and you have other classes here and yeah that's it's a more like it's a advertisement a way to  try to stimulate you for the next sessions until the afternoon and we will really process this  data and produce one like the rival product on top of this lands up data and you can you can use  it for other applications also so inside of the project we are using this data to map where are  the three species for example you can use this data to identify for example and to map the soil  properties and you can there are several applications to build on it but it's important to  emphasize that to depending of the scale depending of the area of your interest if you plan to do  something like that for a whole country for example all this data preparation all this  data immunization it's it's hard working and it's it's hard work and we are providing it as  service for different users and different organizations and possible applications  okay we have more questions or any questions actually  okay okay it's about the difference between the end DVI  why it's subtracting the previous year here so it's the first year yes I think it's related  with yours I can explain yep okay it's more related with the Python syntax okay so  so what I'm reading here is it looks to you like we're subtracting the last year from the first  year right because yeah this the Python syntax is because these this means that so we're not using  the first year we're using every year after the first because the first year has the index zero  and the second one has the index one and then if you say one and then then this the yeah yeah  okay you get it now great so was that it oh yeah yeah indexing from one is is better but yeah  I'm spiteful in Harris from a long line of of working on limited hardware and all of that so  yeah no problem yeah yeah indexes start from one it's weird for me but it's how our manager  think okay so more questions are let's stop and now I think it's time to the lunch okay  thank you all and see afternoon we see each other in afternoon
B1_0_251.png B1_252_3467.png B1_3468_5243.png B1_5244_11411.png B1_11412_12443.png B1_12444_15623.png B1_15624_15839.png B1_15840_16583.png B1_16584_17423.png B1_17424_17603.png B1_17604_17663.png B1_17664_21323.png B1_21324_21527.png B1_21528_21755.png B1_21756_22919.png B1_22920_25139.png B1_25140_25259.png B1_25260_28223.png B1_28224_30503.png B1_30504_31199.png B1_31200_31247.png B1_31248_32183.png B1_32184_32855.png B1_32856_32963.png B1_32964_33143.png B1_33144_36623.png B1_36624_40955.png B1_40956_41399.png B1_41400_42095.png B1_42096_42803.png B1_42804_44327.png B1_44328_44915.png B1_44916_44975.png B1_44976_45047.png B1_45048_45095.png B1_45096_46055.png B1_46056_46931.png B1_46932_47711.png B1_47712_48383.png B1_48384_48599.png B1_48600_50195.png B1_50196_51239.png B1_51240_51767.png B1_51768_53207.png B1_53208_53699.png B1_53700_53879.png B1_53880_55931.png B1_55932_56087.png B1_56088_58715.png B1_58716_58763.png B1_58764_60263.png B1_60264_60311.png B1_60312_60359.png B1_60360_60671.png B1_60672_62891.png B1_62892_63071.png B1_63072_63227.png B1_63228_63287.png B1_63288_63455.png B1_63456_63551.png B1_63552_64055.png B1_64056_66011.png B1_66012_66083.png B1_66084_66299.png B1_66300_66443.png B1_66444_66491.png B1_66492_66563.png B1_66564_67235.png B1_67236_67355.png B1_67356_67463.png B1_67464_67655.png B1_67656_68099.png B1_68100_68303.png B1_68304_68999.png B1_69000_69563.png B1_69564_69971.png B1_69972_70403.png B1_70404_70463.png B1_70464_73691.png B1_73692_75143.png B1_75144_75251.png B1_75252_75479.png B1_75480_77783.png B1_77784_78179.png B1_78180_78455.png B1_78456_78671.png B1_78672_79619.png B1_79620_79703.png B1_79704_80303.png B1_80304_82247.png B1_82248_83267.png B1_83268_84047.png B1_84048_84347.png B1_84348_84455.png B1_84456_84743.png B1_84744_85427.png B1_85428_85739.png B1_85740_88031.png B1_88032_91427.png B1_91428_91571.png B1_91572_92003.png B1_92004_92423.png B1_92424_92519.png B1_92520_92891.png B1_92892_92939.png B1_92940_96011.png B1_96012_99395.png B1_99396_99575.png B1_99576_99803.png B1_99804_104279.png B1_104280_105539.png B1_105540_113243.png B1_113244_114995.png B1_114996_115043.png B1_115044_123587.png B1_123588_126251.png B1_126252_127151.png B1_127152_132743.png B1_132744_133140.png
So, hello everyone, my name is Lukka, I am doing the training session about introduction to spatial and spatial temporal data in Python. So the outline of this workshop first, I would like to start with some spatial referencing basics because I found that a lot of experts that are proficient in GS actually sometimes forget the basics of coordinates, reference systems and what they mean for the data and even myself with actually background in Geodesy and it reminds myself of this sometimes. Next, we will continue with Raster IO and manipulation and after that with vector IO and manipulation in limited capacity because the rest of the workshops are mostly focused around points so we will limit ourselves to that. We will then compute a basic time series and go through some, if there is time, go through some new map convenience and performance utilities. So to continue, so what is actually spatial data? It's data in a spatial context. So when we say spatial context, we usually mean the data is placed somewhere on the surface of the earth but spatial context is more than that because we build multiple layers of abstraction around the earth surface. So we have like first there is a surface of the earth, then you have like a terrestrial reference system, then you have an ellipsoid and then after all that you have like a surface that you project ellipsoidal coordinates on. So yeah, we need all of that to put data in a spatial context. So Rasters are positioned basically with multiple parameters. You have the coordinates of the upper left angle of the Raster. You have the size of the pixels which is actually defined by four parameters, two of which are usually zero. If you see the BND, the BND, sorry. And if maybe some of you recognize this image below the left side, it's from the Wikipedia page on Esri Worldfiles which is basically what's positioning Raster in the spatial context is it's an FN transformation. And vectors are positioned in bits of a simpler way which contains more data but the data itself is written into the vectors. It's just coordinates for which vertex. So yeah, the X and Y coordinates in the context of the geospatial data are a reference to the projection surface. This is also true even for 3D data because elevation data is think of its own basically, usually. So the quality of a coordinate reference system really depends on what you are using it for. So basically for mapping, topographical mapping or web mapping, like navigation maps and that, you would use a conformal projection because it preserves angles. So things retain the shape of how you know them. But for quantifying surface properties, for example, you would use equal area projections like we use for open data science because that's what we're doing here. And why does any of this matter? Because we use data from all kinds of sources and we forget that it might not be compatible with one another. So you can have like, well, if you have continuous data, then if you project it in another CRS, then everything is probably fine because you interpolated, but if you have like, have a nominal data like, for example, land cover, you cannot interpolate it. So you do like nearest interpolation and you have pixels that are largely different shape and if you compare them with another, for example, nominal data, data set, your results would not really hold water. And also with vector data, there's things like old coordinate reference systems, like for example, in Croatia, we have before the current system, we had inherited the system that was from the Austro-Hungarian Empire, which had a local ellipsoid and people are still having trouble to this day converting data from before to our new system because you basically get vectors that are 400 meters away from where they should be and stuff like that. So yeah, all of this being said, you won't have to deal with this during the course of these training sessions because all of the data is in the same system and all of the pixels are aligned in exactly the same way. So yeah. So yeah, let's continue with the practical part. Okay so before we start, if you joined the previous session with Landra, you don't have to do this but if you are just joining us now, please update your notebooks like this. So you go to this directory and just type it pull. So okay, we can begin. So yes, again this session is very introductory. So if any of you are joining this just because you are joining the entire Python track, then sorry if you are not going to see anything useful but yeah, so let's begin. So yeah, we'll cover the basics of accessing basically and manipulating data in Python using the standard tools which are REST, RIO, Numpy, Geopandas and also UMAP which has some performance and convenience oriented functions like Landra already mentioned. So yeah, okay. Let's clear all the content. So let's begin. We will gather some files. I'm not sure if, yeah, Landra already mentioned all of this. So this is, we will use this style, this style ID, this is the data there. We can just go ahead and run this. We'll get a bunch of file paths. Then we can open the second file in the list and then we get the dataset reader. Data set reader is a REST, RIO object that enables us of course to read files and contains file metadata. So if you run the next cell, you see there's a driver, Geotiff. There's the number of bands from the file, the shape of the REST, the coordinate reference system and the transform. No dates of any as well. You can actually see all of this in the REST profile. Maybe don't print it. See the display. No, it's the same, never mind. So now if you read the data, you can put the number, the band index here. Here there's only one band, so it doesn't matter. If you do not put the band index there, it reads all the bands and just if there's only one band, it just adds one more dimension to the array. So we're going to go ahead and do this. So when you read the data, we get an umpire array. And we can do with it whatever we could with other umpire arrays. So we can do this. We can compute with the methods of the array. We can also do some stuff which we'll get to later. We will do that, compute some statistics and then plot with the plotter from your map. So yeah, this is the REST there. I think maybe the picture is a bit big. So if anyone is wondering why the fixed size is only one number, if you're familiar with Matplot Libits because the other axis is computed automatically to preserve aspect ratio of the RESTers. So yeah, here we have a, what was it? It was a blue band, I think. Yeah, blue band. Well, let's make it blue. Okay, so here's the blue band. So now we get to do some stuff with NumPy because NumPy does, enables you to do computation basically parallel under the hood. You don't have to think about parallelization and does it on its own with BLAST libraries. And one of those things we can index, we can compare arrays to a number in a vectorized manner and then index them with what we get from there. So yeah, here we'll compute the 20th and 80th percentile. If you're concerned with the usability of this, don't be, this is not about doing something particularly useful. It's just to go through some things one can do with RESTers, with RESTorio and NumPy. So yeah, we'll compute the percentiles. We'll compare the data with the high and the low percentile. Then we will combine the index with the bitwise or, so bitwise means you compare basically the indices element by element. So you have two arrays of Boolean data and you compare them element by element so you get a combined index. And then we just, yeah, we can calculate the percent of data outside the bounds by, yeah, this code. So yeah. So this is what we get. The low percentile, oh sorry, this should say 20, not 5. And this should be 80. So the low and high percentiles. And this index is, like I said, an array of Boolean values. So we can then use that array. So we will first create a copy of the data so we don't soil our original array and we will index it with this. And say that all of these pixels that are true here should be no data and then plot the results. Now again, reduce the size a little bit. Yeah, okay. So you see now that a lot of the raster is now no data. We can do something that's usually more useful and actually clip the values to the low and high percentiles like this. So we index first with the low index and say that the values should be low and then with the high index and that they should be high and we can then plot the results again. Yeah, again, let's make this blue. Yeah. Okay, actually let's not make it blue because you can see the difference here. I'm going to go back to the previous image and actually just let's plot them next to each other. Data and new data. Let's add another title here. We don't need the F there. Yeah, okay, you can see the title very well. Let's do this. Yeah, okay, so you can see the reduce of variability in the data because the right one is clipped. So there's something we can notice here as well is that this corner here and this cut out here that were previously no data are now actually clipped because no data of this raster is zero if I remember correctly, which is lower than the low percentile and we did not do anything to account for that. So what we can do is this. Restrao enables us to, we could also compare just the data to the no data value and then produce an index that way, but the more efficient way is to use the read masks from Restrao. So we read masks from the same band. We do this as type of to make it a bullion area because you get you in date, I think. If you don't do that, then we can do this again. So we get only the pixels where there's data. So we can compute only on them. So we get data only. We do again, we do the indices. We compare the data pixels to high and to low. We again clip them. We make another new data array. We fill it all of it with no data and then we can fill only the pixels that under the data mask with what is contained in this, the data only array. And again, let's actually let's plot the old data here as well again. No, sorry. No, wait. So if you can see there, no data is now accounted for. So we have made a new raster and let's now go ahead and write it into a file. So first we will create a directory with this part here. Then define the path of the output file. And now we can use Rasterio dataset writer as a context manager as one would with a regular open from the Python standard library, but aside from the output path and the mode, which is the right mode here, we need the raster profile. So if you remember the dictionary that contains all the raster metadata, well, we can just pass that here as keyboard arguments because we haven't changed basically anything other than the data itself. So the data type is the same, the CRS is the same, the raster size is the same, everything is the same except the data. So our new raster basically has the identical profile to the new one. So let's go ahead and write it. You can see we have here dataset writer object and it wrote to the file. So just to check we can use plotters with the file path to see what was written. Again we'd use the size of the plot. So there we have it. So now to some vector data. For that we will use Geopandas library that basically it's Pandas, but for spatial data it has Geopandas data frame basically differs from the Pandas data frame by having a non-optional geometry column on which you can compute geometry operations. So we will go ahead and read the land cover samples that Leandro mentioned in the session before and print some data points and then look at the points themselves. If you ignore the runtime warning there. So the CRS is the same as the raster data. Here we have the number of points and these are our points. So we will not be using any of the data in this data frame actually. In this session we will only use the geometry. So yeah it was mentioned that IPa leaflet is useful for visualizing data in spatial data in Jupyter but if you need a quick and dirty way to compare the extents of datasets this is one. So basically you can import Shapely. So box and multipolygon from Shapely geometry. Shapely is a library that Geopandas uses under the hood but used to now it's used as either that or PyGeos but yeah we won't get into that here. So you create a box from raster extents and that creates a polygon. You then do the same with the points. So you do this. Now the GascadedUnion is if you have a relatively small dataset you can make the entire data frame into a single Shapely geometry so that it visualizes easily in Jupyter and also you can compute bounds on the whole thing basically. So yeah we can do this. We have two polygons and you can see our raster is here in this little corner. And these are the points. The other points are just here and here but yeah we will just, what we will do next is take only the points that intersect with the raster. So we can do this. We create an index by using the points intersects method. So you just pass it to another geometry and you get a Boolean array that contains TrollForce depending on where the geometries intersect with what you have here. So yeah let's print the number of points we get now and let's also visualize the points subset. So yeah these are our points now. Here we can also write vector datasets. So right now we have to do this because there's bug with the current version of Fiona so if you don't do this at least we, I think with our point dataset produces an error so basically we have to do this. So what is important here other than the error wrapper is you just do this points subset to file and then you define an out path and then say which driver. You only have to define the driver if you're not writing an esri shape file because that is the default. So yeah we now wrote that to a file. So onto computing raster time series. So yeah like really Andrew already mentioned there's a helper in UMAP that enables you to read a bunch of rasters multithreadedly into a single array so you don't have to go through all the files manually and then stack the arrays and all of that. So we will go ahead and take all of the spring data set, all of the spring rasters for red and for near infrared for this style that people we were already working on. Yeah here are the shapes. As you can see these are stacked like a multi band image instead of being like this axis being on the front where you have like a list of raster datasets. So yeah we can basically just do this. Raster datasets are arrays you can compute with them directly and this is how you compute an NDI time series from the time series of red and near infrared. And we can also index one point and one pixel well and plot series for it. So yeah that's that. So we can also plot the entire series with plot rasters and I'm not sure if we have to do this anymore. I think Leandro maybe patched this last night but I'm gonna go ahead and do it anyway. Just move the axis that contains the bands to the front. So that the shape of the array looks basically like this. And yeah we are going to plot the rasters now. I think I can increase this a little bit. So yeah there we have it. We didn't plot the entire series because of the 20 rasters. We can try it but yeah. So we can now compute with the series. We can for example compute the difference of each year to the previous one. Which is basically just differencing the series from itself with the shifted index. So we take here the last 19 images and then we subtract the first 19 from them. And then we can also plot the results. So yeah now you can see differences from each year to the previous. Or you could say you want to compute differences from the first year. To this. But yeah. Let's get back to the previous year differences here. And we can use the do-map-helper function to save the rasters in parallel. So yeah we define the file names here. And we write the files. Also yeah we need to put here a path for reference raster. Which then saves rasters takes the profile from if you remember the raster profile you need when writing rasters. So instead of putting that there you just give it a reference raster to take the profile from. So yeah basically we just we unstack the data again to be a multi-band image because that's what save rasters takes. And we also should define no data here because the default is zero and the difference of Ndvi from one year to the next can definitely be zero so we should change that. We also change the data type because the GeoTIF driver does not accept float 16 which is what our data is because we computed from the start with 8-bit integer values. So it tries to keep the memory usage low by using the lowest possible float. So yeah let's save the files. And here save rasters just outputs the file names that we saved data to. So what we can do now with these files we can use the overlay helpers from map mapper to overlay all of the rasters at once with the geometries. So this optimizes things a bit and so let's run it. So is it done? It's done I think. Yeah it's done. So basically give this class a bunch of geometries as a Geopandas data frame and just a list of the files that you want intersected and it does the overlay for you. So now we have time series for each of the points. So you can see the column names are just named after the years. So let's plot the series from a single point. For a single point we are going to sort the columns here because they are not currently sourced by year. And we are going to take a series from a single point and just plot it. So yeah as you can see series at this particular point. So now onto some convenience functions from UMAP. So yeah RasterIO allows for window reading of rasters. So if you do this, if you define a window you can read raster with that window. So only the data from the window will be read. So if you do this, so you need a column and row offset and width and height in pixels. So yeah now we have a 5x5 pixel data set. But what matters here for large files is that RasterIO has to read the entire block of the raster before it reads that particular window. So if a window intersects with multiple blocks then it has to read all of the blocks before reading that window and just outputting the data. So if we look here at the number of windows that our test raster has we can just do block windows and just expand it into list because it's generator. So yeah 125 blocks. So yeah but what we can do with UMAP blocks basically we can read a larger raster by just reading it block by block window and we can not worry about that let it run in the background and just define something that we would like to do with the data from each block when it's read. So we can do this by import mapping from Shapely and take the URLs for a landsat read and near infrared for 2018. This is summer. And so we will define the file URLs and define the geometry. The geometry has to be given here in the JSON schema. So we do that with Shapely geometry mapping from our Shapely geometry. So this is the raster that we, the size of the raster that we are going to read now as opposed to 1000 by 1000 pixels that we were reading before. And this is our geometry. So we import from UMAP Parallel blocks we import the block reader, the block writer and that's basically it. All we have to do now is define NVI as a function that will be computed for each block. So we initialize the reader with the reference file and here all of the files that are read have to have the same transform and the blocks have to be identical. So this is for reading those kinds of data sets. So we give it the radio as the reference file. Then we initialize the block writer with the reader. You can also omit this and you can just initialize the writer and then that will initialize the reader with the first file. But just to illustrate because the reader is then reusable, it has to, it initializes by computing geometries for each block. So it can take a while so if you are going to reuse the reader class, the reader object rather you should initialize it on its own. So if you run that and wait a bit. So check. Oh, it's not done yet. Okay. I should run this. Yeah. Okay so if it's done. So we did, we defined our NVI function here. We defined an out file. Then we called writer write with the source paths and the destination path. And also the geometry, the function that we want to compute and should give, of course, the no data value because again I think the default is zero. The d type and now let's see if our files have written our file. Sorry. Yeah, so this is the same block but computed from very large files that are hosted in S3 buckets. So you do not have to download the entire file to work with it basically. Yeah, so the final part here, I would like to talk about some convenience abstractions around the catalog of the data produced and also used by Open Data Science. So as you know as Leander has mentioned, we have the catalog in the geo network but from a Python standpoint you would have to use CSW from OWSlib to access that directly. So here we have a catalog object that we import from EOMAP data sets. And then we can search with it. So here the catalog also has a fallback copy of the basically of the entire catalog as a CSW that's distributed with the code. So if basically if maybe you're having connection troubles and the unit per can also be a bit slow you can just use this and that's what we're going to do because it's faster and the fallback CSW is updated so there shouldn't be any difference between the two right now. So yeah, here we have a resource collection. So the output is grouped by metadata so all of the resources that have the same abstract theme and title so this is NDIA Landsat Quarterly. They're grouped basically but you can use the results object as you would list so we can just do something like this basically results. So to get the first three of them. So yeah, these are and if you look at only one of these. This is so these results objects are basically strings so you can use them as URLs for files to use with Rasterio and the UMAP helpers but they do contain additional metadata so if you want only the metadata objects you can do this. So here is the metadata. But yeah and also there is the metadata of the collection object. If you want to access the group, the metadata groups of results so if we have zero that gives us the quarterly Landsat NDVI. We index the one. We get this and oh no there's no three. Okay, this output is a bit long so we'll scroll down here. We could also use, you could exclude the results with the search so we don't have to sift through them manually later on we can just exclude with the exclude keyword. So if we exclude the trend we just basically get this. So yeah let's check results. So yeah, there's basically only the first one. Again with the long output. And we can also search by year. So here we get the quarterly Landsat for 2019. If we want to say search by multiple years so we can here input range for example from 2015 to 2020. This is non-inclusive so we should be okay. Then you have, then year is added to the metadata and you get groups of resources by year. So if we do this. So that's now just the first group. So only the rest is from 2015. So yeah we can use now these resource collections to read directly from the files. So we'll use the read rasters. And read rasters also takes a window object so we will define a window like this. So yeah the rest of the window sub module has this from bounds function which takes basically the bounds. So we will give it only the bounds that we were using so far so the extent of the test raster that's included in the VM. We will, we also have to define a transform. We can also just read from one of the files so we get the first file. We just do this. If you remember it's what the rest of the data set reader contains of all of the metadata like this. So yeah you define a window here. You just basically use read rasters as we did previously but you just say you use the spatial bin keyword argument to define a window. And this part here is basically the catalog result resource should behave like regular row strings but there's currently bug with that so we just convert all of them to strings. So we run this. And then we can plot the results. Okay this is a bit too small. Okay let's instead just plot. Well we can do that here actually. Oh sorry. Yeah so these are the windows read directly from the rasters that are hosted in S3 buckets. So okay so this is the end of the plant section so we have some time left if you have any questions if you would like to suggest that we try anything or. Are there any questions on zoom? No? Okay. Andra would you like to use the time left to go over the things from your last session that you had to skip over? I can just explain it. Okay. I think that's a good idea. Okay yeah. And yeah so you already covered some parts but I can open and explain a bit. Okay great. Thank you. So just a short recap. So we presented like how basically to access the data, read and write parts of it using the local data that we have in the virtual box like in the machine. And some of the remote data also. So this is actually the first part so to really read the points, read the raster data and in the afternoon we will use these functions but now in the context of spatial and spatial temporal machine learn. So all these functions that we present it's mainly to implement these input and output operations like deal with these GOT files. And I would like just to in my last, in my presentation the first part of this morning. So here you can see that I didn't present it but it's the same method that we use here to read these three images. You can use the same approach. You just need to change which are the file system, the file names and with these file names, these function files, it's like a utilitary function of EU map. It will sort all the files in the chronological way because the file name, you have the date inside of it. So mainly this function organize all the files in a chronological way and you can read all these files using this red raster and it's fully optimized so it reads in parallel. You can control how much threads you can use to read it and you read here all the 84 images. So here for the NIR like the infrared band of the Landsat so this is a totally harmonized product so you have it for the whole Europe and you can read it and here you have like a complete time series with four observations per each year. So we have 20 years and here it's actually 21 because it's 2000 to 2020 multiplied by four and what it's interesting so here it's just like a noon pile array with the last dimension it's the time. So each part of this array so it's like a data cube with each part of this last dimension one date but with this function here you can see I'm just plotting like a time series and I'm getting like a random pixel here and you can see the time series value so how for this specific pixel this pixel it's one position in this data cube structure you can see how it is changing. So tomorrow in the cloud optimized geative training session we will do it but with real data so we can actually access this temporal response for NDVI data directly through the cloud so it's here again you are reading the data it's just local data it's just part of Europe because we prepared it specifically for the training sessions but you can use the same approach with the same methods and you can point this to like a real data that we are hosting in the cloud and you can plot all this time series for NDVI in our app. Example and here you can compare so this is actually the structure of the data so you can compare it across different years so this method here it's just like a function to help plot the data using like the plotter the plot haster that we saw here and you can send like two years and you can compare how the image looks like so the four seasons across different years and here it's just like a operation between this this image so it's a sort of reduction so considering all this data cube you can do multiple operations and in this in this case it's just like a median you can use like maximum calculation minimum so basically it's a type of a reduction in this data structure so you have like two dimensions for the space and 84 dimensions for time each of these dimensions are specific for one season and for one year and you can derive multiple statistics about it so for example here it's you this code it actually calculates like the median value across all these images across all this time series you could have like the maximum value the minimum value so it's a easy way with the work with this data and here it's just an optimized way to access it but all the functions available in umpai you have available also to process and to deal with this data because we are creating this EU map function just to provide some helpers and to optimize for example specific parts of workflow but it's fully compatible with what we already have in python numpy scikit learn as you will see like this this afternoon so yeah this is the last part of the presentation and you can see it's possible also save the raster files and you can open it in the quantias and it's a nice way because when you save a file you need to you are reading a data put it in memory so you have a numpy array you can process and do multiple operations but to save it as image again you need to set several parameters what is the projection system and what will be the like the the data format and things like that so this method save raster it's actually used as profile one like base image and copy all these geo spatial parameters for all these images separately and you can save multiple files like that also so yeah that's that's it's mainly what we want to present this morning it's more about what data we are processing how we are reading and how we are writing and how you can access it using python so it's a preparation for the afternoon because in the afternoon we will use this data to do a proper like machine learning model and classify it and this landsat data primarily we will classify it in different land cover class so yeah I don't know if you have some questions so otherwise we can stop the session have a lunch and prepare for the afternoon sorry okay okay so let's okay let's go back here and I think if we don't if we don't have like more questions I think we could start with yeah I can actually I can execute it and because and we can yeah okay so I will open mine so I was rushing but as we have time now I will I would like to explain about this data set preparation so it's a it's a nice context for the machine learning model so mainly what we are this is mainly the data that we work until now so the landsat and R&D images so you can see here this is a nice publication about the data set but mainly we are downloading it from University of Maryland as I said and we we are getting all the like this 16 day style so for the whole world is this product provides like one integrated image every 16 days and when I say integrated image it's put together all the three landsat satellites that are available like landsat 5, 7 and 8 so it's a nice product it's a harmonized product which makes all these pixel values and when you think that there are three satellites different getting data from the earth surface of course you have different problems and it's really challenge it's really a challenge put all this data together so this is a harmonized product it's it's nice but even considering it you still have miss data and you still have like clouds so what we are doing we are getting this this data and we are organizing it for each season and mainly we have the date intervals here so for example between March 21 of March and 24 of June we will have like several images maybe six seven it depends in winter of course we will have less images because you have more problems with cloud and but it it changes across all these seasons but for this time frame for these dates you will you will have like multiple images mainly what we are doing here is we remove all the clouds so there is like a mask there is an indication of which pixels are cloud which pixels are not cloud and we are using this mask to really remove and clean the data but of course if we remove the cloud so this satellite it's not it's an optical satellite so it doesn't see across the clouds so you will have a value and a gap there so and if you have a gap if you have a miss data and no data it is a problem to multiple applications so it's it's it's better it's actually it's mandatory put some value there to use this image in this format but if you think that okay maybe for this interval I I will have like then this data will have gaps so there are two options here you can expand this window so imagine it if if I have like value pixel values here cloud here I can expand my window so I can get all the images available I don't know in 2020 for example and maybe some of these images will have like pixel value so it is not guaranteed that I will have like observations because it's a combination of when the satellite pass and if we have clouds or not what are the climate conditions but you can expand it but of course if you expand it you will have like a final image with different with pixels derived from different observations and it could mess the data so it's it's why really it's important harmonize all these observations so and but in the context of the project we decided to use we decided to use other approach so mainly we kept we kept this seasons and we get feel with you we get fill it using like different years and neighborhood seasons so and I didn't prepare the tutorial about it but we have available in the official documentation I will just show to you but we in the EU map documentation and in the official documentation we have these tutorials so and here it's like a it's exercise of how we are implementing the gap filling so this is actually the code that we run so and you can see here where is the original data and where we get filled the observations so mainly it's it use like a moving window to get observations from different seasons and from different years and fill the gaps of course you need to be aware that you are introducing like a new observation there that and in some applications it could introduce some problems so imagine it like in a case of a deforestation if you put like a gap filling maybe you are not see a real deforestation because you are using like a value from the past for example but at some point we will see that and because at some point the satellite will pass through this part of the earth to that part of the earth and you can you and you will have like a real observation for surface and you will see it so it depends of the application but here as we are mainly we are dealing with data to the past so past data since 2000 we have like some observations up forward also so what we are doing we are just using like the medium value from past observations and from future observations so for example to get filled I imagine 2020 2010 we will use observations from 2012 and from 2008 of course we will have limitations in the end of the year in the end of the time series like 2021 so and we are preparing like a publication to explain it but for the context of the project and to produce like land cover maps to the past this data it's it's it's prepared to do it and of course we have four seasons and maybe one of these seasons has actually observations without gap in the most part of the case that's that is the case actually so maybe we are gap filling for example the winter image because there are too many clouds in the season and but the summer image we have observations so it and we are using all the seasons inside of one year to produce the land cover maps so this problem will be minimized and and here you have other approach we did test three of it and one more that it's not implemented here and here you can see this is more like a linear regression between all the available observations in the cross time but you can see that our approach here it captures better the decisionality and we see these spikes because we actually we just have four observations per each year but if you think that this data it's a 30 meter data so and it was prepared for the whole Europe so we removed all the clouds and it's a nice time series and it allows to develop several applications on top of it and so let's go back and we will use a lot this data and so just to explain like the last part we did it for three percentiles so mainly inside of each season imagine you could have like four observations three observations maybe you have just one observation that's it's the case for the winter in the most part of Europe but we are calculating these three percentiles for all the seasons so it's a it's this approach it's kind of it's a straightforward way to produce to organize these different data availability in three ready to use images so for example if you have like just one pixel and for sure this three percentiles will be the same but if you have like on the Italy for example you have a better observations and a better data availability you have three different observations for each these three percentiles because you will combine all these pixels and you will produce it so the output here and of course it depends of the data availability but it's a it's a approach that allow us to generate this data and to prepare it for the machine learn and the assumption here is when we train a machine learn we will use points from the whole Europe so we will have all these different situations so we will have points with three values constant values in the three percentiles inside one season and we will have values with like multiple person with multiple values for it is different values for these three percentiles and the machine learn the model will manage it so that was the assumption here but again there are multiple ways to deal with that you can we could use like just one image for the entire year that that's a way I implemented it in the past but here we really preferred to generate the seasons and get through it and and it's it was more like a methodological decision but of course there are different strategies that could be implemented the important part and I I'm repeating it because it's not it's not these satellite images they are not like if you think in this Landsat data in the Sentinel data as they are available now they are not ready to use you need to do a lot of processing and you need to download it you need to remove the clouds you need to do some a lot of harmonization to use for example different observations from satellites so and there are there are like specific projects to do it so here we are using one of these products that already implemented it and we are building like a new kind of additional product on top of it and if you use directly the satellite images so you have a lot of work and probably during the machine learn modeling you will have you will have problems and because this even for example from the same satellite you can have like different interference different like noise from the atmosphere effects and things like that so if you are really focusing develop like models and use this data it's important to have a clean data harmonized and ready to use otherwise you need to implement that and that's the reason of why we are doing it and we are putting all this data available so to allow other people other researchers and other organizations to develop applications on top of it and of course we are improving this product so now we are generating a new version a version 1.1 and it's still this product it's not it's not perfect you can see like some artifacts still but we are really improving it and during the year so and we have the capacity to rerun it so I will as I'm talking about it I will open like one RGB image here so we can use the data portal to see it on top of my tutorial it's the easiest way to check these composites and if we go here so this is a public URL you can access and this is the RGB data and I will use the zero capacity and I will zoom in some place you can use this to yeah we can try to find one and so yeah this is nice so and we are here and you can see this is this is not like a Google or image because of course each pixel here has 30 by 30 meters so if you zoom in it it will not shift for a better image as we are as Google Earth do you have the Google that not here is the big map but it's just to show the difference but this image it's really enough to do like several environmental applications and so but what I want to show through here it's this is the summer image so it's the best mosaic that we have available because of the data availability and this is the medium value and here you are seeing actually three bands like the red green and blue and if you go back here you can see like from all the years you can see some changes and you can see still some artifacts so it's there is like rest of like some noise some artifacts because of course this processing it's not perfect we are trying to minimize the number of artifacts but sometimes there's problems in that there are problems in the cloud mask so it's really difficult remove all these problems but what it's important is we are really like creating an optimized workflow to do it and to update it and for example if we find a better way to identify the clouds we can run it and we can process again and using all these infrastructure that we are presenting here so what it's nice so we are here in Farnigin and this is Fenendal and if you see all these images you can clearly see some of the urbinary expansion here and this is only possible to see because we can have we have one image per year here and for each year we just have we have four images actually here it's just the summer but for Google Earth you have this for some areas you have multiple images but it's not really you can't do analysis with that so with this data as we are we presented today you can load the data in the memory and you can process it you can derive some statistics you can calculate like NDVI index index and for the NDVI we put it all and now I will turn it so if you look here for the NDVI you have the four seasons and I will turn on the opacity and you can do the same thing and see how the NDVI changed and as we know this now this area expanded so we had an urban expansion here and the NDVI will change because in the 2000 and the early years probably it was some kind of crop or pastry area and later it was converted to urban area so the vegetation index changed across it and if you you can check the values here so you can clearly see how it's dropping it and so probably in this specific in this specifically in this pixel maybe the conversion occurred here so and again this is a data that we prepared and we can try to identify it here yeah so you can you can see that you want considering one year to other year I will stop it here so yeah you had an expansion so we use it this data and we produce it as I presented but I will show it because it's a nice example so we produce it this land cover maps so we prepared the Landsat data we harmonized we remove at the clouds we put it available and we calculated this NDVI index and on using this data so it's not really raw data it's an harmonized product gap filled and ready to use but it's just data right it's just like we just have a spectral value like a pixel value that represents some some response of earth surface but it's not there is no like directly mean for this pixel I don't know if it's a urban area I don't know for example if it's a crop if it's a pasture and use this data we actually produce for all these years we produce it like a land cover map and we actually classify it using these spectral values we classify it as my internet here it's not helping you see yeah now I have it so here it's like this is the result of the machine learn model it's so you can see that this area were like it was a pastry I think and so what it's interesting here if you look the image you can identify it but yeah this because we all most like the most part knows like what is a urban area we can see that maybe the we can see the building growing up across the landscape but for a machine learn model and to produce like a product on top of it we really need to train a model we need to provide training samples so there is a whole workflow on top of this data and if you put here so here you can have like you can see all the land cover class and when it was converted to the urban area so and this this is like the result of the land cover classification we will do it on the afternoon so we will produce something like this as output for a specific for the tiles that we we have and using the same the same workflow and you have other classes here and yeah that's it's a more like it's a advertisement a way to try to stimulate you for the next sessions until the afternoon and we will really process this data and produce one like the rival product on top of this lands up data and you can you can use it for other applications also so inside of the project we are using this data to map where are the three species for example you can use this data to identify for example and to map the soil properties and you can there are several applications to build on it but it's important to emphasize that to depending of the scale depending of the area of your interest if you plan to do something like that for a whole country for example all this data preparation all this data immunization it's it's hard working and it's it's hard work and we are providing it as service for different users and different organizations and possible applications okay we have more questions or any questions actually okay okay it's about the difference between the end DVI why it's subtracting the previous year here so it's the first year yes I think it's related with yours I can explain yep okay it's more related with the Python syntax okay so so what I'm reading here is it looks to you like we're subtracting the last year from the first year right because yeah this the Python syntax is because these this means that so we're not using the first year we're using every year after the first because the first year has the index zero and the second one has the index one and then if you say one and then then this the yeah yeah okay you get it now great so was that it oh yeah yeah indexing from one is is better but yeah I'm spiteful in Harris from a long line of of working on limited hardware and all of that so yeah no problem yeah yeah indexes start from one it's weird for me but it's how our manager think okay so more questions are let's stop and now I think it's time to the lunch okay thank you all and see afternoon we see each other in afternoon
to start with some spatial referencing basics yeah, we'll cover the basics of accessing also visualize the points subset. So yeah these of the data produced and also used by Open Data could exclude the results with the search so we to access the data, read and write parts of it to really read the points, read the raster data compare how the image looks like so the four