Handling and analysing vector and raster data cubes with R
CC Attribution 3.0 Germany:
You are free to use, adapt and copy, distribute and transmit the work or content in adapted or unchanged form for any legal purpose as long as the work is attributed to the author in the manner specified by the author or licensor.
Pebesma, Edzer
Vector and raster data cubes include vector and raster data as special cases, but extend this to vector time series, OD matrices, multi-band raster data, multi-band raster time series, multi-attribute vector or raster time series, and more general to array data where one ore more dimensions are associated with space and/or with time. Examples come from pretty much all areas dealing with spatiotemporal data. This tutorial will go through a large number of examples to illustrate this idea, mostly focusing on the packages stars and sf and those supporting their classes (like tmap, mapview, gstat, ggplot2).
 Anyway, in this GitHub page, you will find the Rmarkdown file. And the Rmarkdown file, you can download, and then you could go to download it somewhere. And then you could go to the place where you downloaded, download the whole repository, like I have it here. Oops, a lot of things there. But there's also the R&D file. And you can open that with RStudio, and then you could theoretically go through the code sections and basically run the whole thing. So you can start there, and then whenever there's a code section, you can run this, and then go to the next code session and so on. And run all these things. And hopefully this command might also work. This command actually downloads a somewhat larger data package with some real remote sensing image in it. So that's a package that has a size of about one gigabyte that you can throw away afterwards. But that I use because Chrome packages can only have a certain size, right? So you're not supposed to send your one gigabyte packages to Chrome and have it distributed there. So I set my own repository where you can install things. I think it also works with HTTPS if you prefer. Usually you can do it just now. Although Tom says not do it, but it typically it works, but just goes a little bit slower if multiple people do it at the same time. So that is where the bigger examples come from, but I will go through a lot of smaller examples primarily. First, we'll say something. So what are data cubes? So we looked at a little bit that this morning and there are essentially rectangular two-way tables or multi-way tables where there's not a structural difference in the data. So we think of data frames also as tables. The data frames are essentially collections of records, right? Like a data frame table. The way they are stored is that they're collections of columns, right? Because columns are vectors. Vectors can be of different types. So you can have a vector of numeric, you can have a vector of date, you can have a vector of logical, you can have a vector of factor, right? The categorical variable and essentially have a collection of these. We think of them as collection of records. The way they're stored is that they're a column store format. But in essence, it's not a two-way layout table with the same thing along rows and columns. It's if there's heterogeneity along the columns and there's homogeneity along the records. And that is basically why tables or the data frames are thought of as one-dimensional structures. So one-dimensional tables are also one-dimensional data. Of course, the sort of the analogy of a cube fails there. But if you think arrays is being generic, having one or more dimensions, it is a, the table is a one-dimensional cube. So typical data cubes is that they have two or more dimensions. Then it becomes interesting. The cube metaphor is of course three-dimensional, but we can also have four or five or eight-dimensional cubes, although you don't often run into them in the wild. So it is pretty hard. Actually, you can find five-dimensional data cubes. You have to do, you have to really make an effort and go to era five, whether forecast data, these kind of things to find them. So tables are a special case of data cubes, but they're essentially one-dimensional data cubes. And raster layers are also thought of, raster sort of raster data are also thought of as data cubes because they have rows and columns and then cell values as the attribute. Then it's the question sort of what are, what are arrays, what are the cubes and what are, what do dimensions of data cubes? Yeah, so I was talking about rows and columns of data, but the typical thing for spatial temporal data cubes is that they have spatial properties and quite often that they might be sort of aligned with some kind of spatial dimension we have, like latitude or longitude, and have some regular, they might be ordered. Not a property, and they might have a regular discretization, in which case we have a regular raster. Right? So, so it goes on and sort of specializes and specializes. The things I will use in these library stars that you can, the best is to install it from get up here straight away, then you have sort of the development version. It's also with crumb. As I mentioned, the vector data cubes, an example is, is a panel data where we have by state and by year we have a number of socioeconomic data. This is a very common way of aggregated data for socioeconomic data, for health data for instance, where you have time series over, over polygons comes off in this way. If you would plot this, then this is said, this is the tiny layout where everything is in columns and basically stacked and you see that, that year cycles fast and state cycles slowest, that doesn't really matter. If you throw that into packages like ggplot, they really like that, right? They like the tiny layout and then they can work with that. You can even plot a raster of this data where you have essentially here, I took the state as the, as the y variable because then you can still read the labels. I think if you put it along the x-axis, you cannot, you cannot read the label and put year along the x-axis because we think often there's sort of time variability in the, in the x direction going to the right. If you plot this data this way using color, so we have a two-dimensional table essentially plotted with color intensities as the, as the value for the, for the public capital variable. And what you see here is of course, it's kind of the most boring thing. It's kind of says, oh, there is only sort of differences between states. Yeah. So this is thing one and then you start thinking, yeah, well, this is public capital. Let's look at which state is, oh, there is New York, which is very high. Oh, there's California that is very high. It's very unhelpful that this thing is, oh yeah, it's ordered from top to bottom to top. But then you see it's kind of the differences between states, which is of course like, okay, yeah, there's more people in those states. So there is also more capital in those states. So it's the trivial signal that you find there. And then you're interested in basically this signal kind of overshows every other state. First signal, yeah. So it might be more. What would be the temporal variability? Are there, is there any change over time? You see, this state goes from a little bit from dark to light, but it's very hard to discard because most of the color range is sort of taken up by the differences in states. So we could do that, and I'm doing that in an array approach where I convert this data set into a stars object with the ST underscore S stars. And this state just takes a table and then it takes it to my default. It takes the first two columns to work out the array geometry. So it says this is the default of this STS stars method for data frames is to date it to first columns. And then y decreasing is false, meaning it doesn't sort of refer to y ordering, which would be the usual case for us today. So I end up with all these variables that started variable three, as you see here, region, private capital and so on. And I move now the first two variables into a dimension table. So this is the dimension table where we see that we have 48 states. We have no offset in delta for this thing. A references, I will talk about what that means. Nothing known about point, but we have values, basically the values for this dimension are the county, sort of the state labels. So the state labels are preserved there. We have another dimension, which is year, which ranges from 1 to 17, which has an offset in delta, meaning it has a regular discretization. So it starts at 1986, and then it has a step size of minus one. The minus one is here, sort of it reverses the order for reasons that you would typically do that by for raster data, which it shouldn't do here. So something went wrong. But anyway, it's fine. But we now have sort of all the things laid out as a two by two table. We could, for instance, work on that table by applying things on a certain dimension. So what I wanted to do here is for every row that I see here is take out mean effect. So because these states differ so much in mean, if I want to look at temporal variability, which is my main interest, I would first take out some kind of state effect. And I do that by dividing by the state mean. You could also subtract the state mean, which is not our approach. I divide it by state mean and say, OK, what I want to do is for every state, so I want to apply for something, a function over a dimension of this array and do that over it first. And then I have a state and then apply functions and then the function that takes the vector that means the time series and divides it by the mean of the time series. And then I retain the two variable names here. And then I get this thing out that just reversed the order of the variables, but still has the same dimensions. And you can see it's now entirely different values, or the values have been normalized by their state mean value. And if you plot that data set and also change the order of x and y and fill it with a public capital, then you see suddenly that all of the states have the same kind of mean level, mean color level, and that we see that certain states sort of go down and certain other states go up. So suddenly we see where we first saw only sort of the main effect of different states, which is essentially the size of the state. We now see that the temporal development by state is also varies pretty much. So it's an entirely different mess. So this is something you can do by saying I want to do something for all the time series, right? So for all the time series, I apply function. This is functional programming and do that over a certain dimension. And that is something you can do with apply, apply works on arrays, and sd apply works on stars arrays. And stars arrays are just arrays with a little bit more semantics. For instance, for instance, it has semantics that we have a regular time series here, although it doesn't know that these are years. It doesn't know sort of how the time information is there. Yes. Yeah, always. So in this plot, the states are different size. And if you would like to plot that also, that you can see the proportional size of the state. That would be possible. Well, you see the proportion of the size of the states. In this plot, by the variable capital. So this is the size of the state expressed of the public capital of the state. So of course, there are states where people are rich and people are poor. That's not in here. This is just the amount, the amount of something, right? And bigger states, this is the size of the state expressed in this. But here state is just labels, right? So we just are interested in the differences between these. Right. And I was here. So this is spatial temporal. It's just that neither time, so we see now that we dropped the year numbers. So neither time nor space is referenced, right? So time was just numbers. We need to know that these numbers refer to years here. And space is state names, right? So we just labels for space, right? And ARC can deal with those kind of things very well. But the thing is that at some stage we want to make plots and then we need to associate these. If we want to make plots of maps, then we need to associate these labels of the state with the geometries, with the polygons of the states and then do that. And then we get real spatial temporal errors, right? So we can do that. And we'll do that with an example of the North Carolina Sitz dataset. So that's a very classic dataset that comes back in a lot of examples. And had a long sort of v yet in the SPDAP package here. And that contains data, epidemiological data. It has several called birth number of births, then the number of certain infant deaths. So this is the disease cases. And then another variable, the non-white births, because there wasn't the hypothesis that the disease might happen more or less with white or non-white persons, and was collected for two years. So we have essentially counties, right? For North Carolina, we have 100 counties. We have three variables. We have two years. So that is like a table with a three-way layout. County, variable, and the time, the year. So we can do that. Here is the table. If you download them and plot them as a table, I sort of used read underscore SF here so that you get automatically get the thing printed as a table. Some people find convenient, which I find convenient, actually. And now I'm going to make this into a matrix. So this is pure base R. So I'm going to select the variables here that I'm interested in, the six variables, because you can see that these variables are spread out over columns. So we have the birth, certain infant deaths, and non-white births for 74, and the year is glued to the variable name. This is not uncommon. We have things repeated in time, and they're essentially different columns in the table because they were associated to the same features. So this is, in that sense, is non-tidy because it's something that's distributed, kind of hidden in the column name. But that is also not very uncommon. So we take these six variables, and then we're going to convert that into a matrix. So that will be a matrix with six columns and with 100 counties. But then we're going to change that in 100 by six into 100 by three by three, just by setting the dimension of this matrix. And now this works because the ordering of the values in the matrix doesn't have to change. So you need to know a little bit about how this works. And what I also do is that I set the dim names. So that means that for each of the values of a dimension, I give labels. So I give the county labels to the variable county, which is the first dimension. I give this birth, certain infant death, known white birth to the second dimension. And I give the year, as you can see here, 1974, 1997, to the third dimension. So then I have a three-dimensional array, 100 by three by two array. And I'm going to cover that into a star's object. And I get this, right? And then it sorts this out. And it says, OK, I have a county dimension, where here the counting names. I have this variable, which says, is this about births and infant deaths, known white births? You could argue this should not be a dimension, but it just threw it into a dimension because of population numbers. And I have two years, right? And I have one variable called population. So I have 100 by three by two. So there's nothing spatial temporarily on this. It is just a regular array with labels. And the labels end up in the dimension values if I throw that, if I convert it into a star's object. So it's the way the dimension names end up is here. But there's still nothing about references. There's nothing about geometries. There's nothing about knowledge that years is actually a time variable. That will all come a little later. What I know now is that from sort of a data cube with labels, with label dimensions, I'm going to convert this into a spatial data cube by setting the geometry of the North Carolina counties to the values of first dimension. So I have this star's object here, which I had here. And I have here the values, which are the county labels. And I'm now replacing that with the actual geometries. So you see now that the values are no longer county labels. They're gone, but they are geometries. They're abbreviated simple feature, multipolygon representations. So the geometries are actually in there now. So now I have an object that is spatial temporal. What you also see is this little thing is that there's a reference system now, which is the coordinate reference system of the coordinates. So these coordinates only make sense if you know what they mean, right, as these degrees latitude longitude. If yes, to which data are they associated there, the data associated to the North American data from 1927? Or are they projected in what is the description of this projection? That all goes into this reference thing, which is the first spatial data, which is the coordinate reference system. So it has spatial data now, and it has a spatial coordinate reference system that ends up here. So this is something that standard arrays cannot do. You cannot sort of go further than just dimension labels. So textuals, it's like a little bit like what I think NumPy also have labeled dimension or x-array and so on. But if things are spatial, then you need to put a little bit more information in there to make this understood and intruparable. And these are kind of ugly things, right? They are polygons, multipolygons, holes and everything. So the whole terrible thing that essentially is solved by SF package is in there. You don't have to worry about that. So what could I do now? Well, I would like to plot them. What shall I plot? I have to think. There are 100 counties. So I could make six maps out of that. What I did is make three maps, but take the sum of the years. So I say I'm going to look over the first two dimensions. I'm going to look over the counties and over the variables and then take the sums of the third. Here's the sum. Here's the function that I call as a third argument. Apply sum to the third. So I mean that I'm summing all the 1974 and 1979 cases and I get that up. That means that I reuse this dimension so that I end up with a 2 by 3. So a 2 by 3 layout. And if I plot that, the counties are plotted for each of the three variables. Whether this makes sense or not, but it's an exploratory plot. And you can still see some colors, which is convenient. You see some variability, although these disease cases are extremely low. They're pretty close to zero in the burst. It's kind of one thing that happens in 100,000 cases, but you still see some variability because the default plotting scale, the color range is cut into quantiles. So it takes essentially the whole cube, all three maps, and then it looks at whatever it takes, like 11 different colors and takes the quantile values. So you still see a bit of variability in the lowest map as well as in the highest map. Although the value, you cannot read sort of what the true color breaks are here. You would have to compute that a little bit or choose a log scale or something else. You could do all kinds of things. So if you want to analyze this data meaningfully, then there is again this problem that we have states with more disease cases and then sort of counties with more disease cases and counties with less, but then also there are counties with more people and counties with less people. So these are just accounts, right? So we are anyway looking, we look at the disease cases. We are looking at a population distribution. So you want to take that out. So you want to sort of have the relative numbers. So we are sort of relative to the amount of births that happen. Where are the sort of the more certain death cases related to that? And the usual way or one usual way to do that is to go from raw accounts to standardized incidence rates, serve values by computing the incidence over the population and dividing that by the mean of that ratio. So the mean of the total population for the whole region. So this is a job that we basically can do over, that we have to do over these three variables. So we have to do that for every year and then basically relate the SID cases to the birth cases. SID is the CIs and birth is the PI's. And that is what I did here. So what I'm going to do now is to move, to split a dimension over an attribute. So like a data frame, we can have, you know, we can have several variables. In this case, we have one array, 100 by 3 by 2 array, which is population. And it's a mix of birth counts and of disease counts. But if we want to sort of do the ratio of the disease over the births, we can as well split that over attributes. And that means that we now not have one population variable, we have three variables. These three arrays are now 100 by 2 each of them. So 100 counties have each of them have two years. So that's another way of looking at the data. There's not a correct way. There is just whatever is convenient. And you can split and you can merge. So you can throw sort of layout dimension over different variables as you would when you create a data frame. This is a list of arrays. So it is kind of a data frame that is smart enough to understand dimensions. And you can split them and you can merge them again. So we split them here. And what I did then is over-disclated data frame, I computed the sum just over this list. So this data structure is just a list. So I can do an s supply, which is a base r function that just looks over a list and see what's there. Average element of a list, I computed sum. So I get the totals, the total birth, total sits. This is the total birth, total sit. So that gives me the denominator here. The totals, sum of ci and sum of b. So I have the sum of birth, sum of sits, total, whole data set. And I can compute that ratio. So that gives me the normalizing constant, so to speak. And then I'm going to use that to compute for every county and every year. Now I want to do this for every county and every year. So I'm going to apply to every county and every year. I'm going to apply, compute the ratio of ci and pi. So this function is an anonymous function. It just gets a length three array, right, with the birth, the sit, and the normal birth values. And I just want the ci and pi, so the sit value over the birth value. So I'm going to just compute that ratio and divide it by my normalizing constant, which was the mean ci or pi for the whole data set. And that gives me the standardized incidence ratio for this data set. So it's a few lines of code, and it looks very complicated, but you just do things. You do it now for every county. Now do things for these variables. Now do things for every year and so on. But keep track of that. If you just, you know, as an exercise, do the same thing for a regular table. You can also do that. It's not rocket science. It is just, in my eyes, it's much more complicated to do it the other way around. And the advantage here is that we end up with two maps. We just don't end up with a table that we then, again, have to link to geometries, but we end up with two maps because the geometries are retained. And the years are retained. And I have CER values for the two years. And those can then be compared. So these two maps are then typically the input for further analysis and do questions like, okay, is there a difference between the two years? Is there some kind of a regional effect at the regions that are higher or regions that are lower? Is there some kind of clustering in this pattern? Is there, yeah, so this is typically the input for spatial statistical and epidemiological modeling with some kind of lattice data models. You see here that I took a dedicated color scale. So I took a dedicated palette, which is a diverging palette, which goes to colors from the center. And I took dedicated color breaks here because the value one is special, because the value one says that I am identical to the mean. So white cells are pretty much identical to the mean. Red cells means that I have incidence values that are higher than the average over all the countries in over the two years. And blue cells means that there's cells that are lower. So the value one is kind of identical to the global incidence ratio. That's a normalizing constant. So this is much more informative than focusing on maps like this. These are the raw counts, which contain, yeah, the SIP map contains the signal of the birth count. So it's very, it's impossible to see whether there's higher incidence rates. And as I said, the advantage of doing it this way is that you can work on dimensions directly and you end up with maps that where the spatial data is still retained. Yeah, so are there any questions related to that? Yeah. So, yeah. You have the June stars. You have the June stars. That doesn't seem to work for you. What do you mean? Oh, that should actually work for this. Oh, good question. Yeah. I think it should work. Well, let me just, yeah, I think it should work because it splits the on stars sort of looks at my roster data and my factory data and then calls the MSF in case of factory data. So it should, but let me sort of you, you're catching me off guard. Obviously. Yeah. There's always a lot of things that don't work that might might work. But yeah. Right. And then can you, you can go on and you can call a map you map you as a currency. It's a little problem that it didn't work with Mark down very well. So I, I, I, I'd comment, but you could look at these plots interactively on a, on a leaf with background and see where these things are. What the names state of the name is, if you click on a particular county. Okay. Let's, let's take a second dive into vector data cubes. So, so yeah, so vector data cubes are very convenient for cases where you have not just sort of vector data with, with attributes, but kind of structuring the attributes in the sense that you might have a time series or you might have a set of variables that is in some way related and you want to work with that on that as a, as a dimension. Let's look at a more complicated or more challenging case where we have more than one space dimension and the case there is the case that I worked out here is that of origin that destination matrices and that case we worked out here in 432 of, of a book that is slow, very, very slowly progressing that builds on data set that Robin Lovelace, one of the authors of the geocomputation with our book shared. And that is an origin destination cube from the Bristol region. And the, the, how am I just going to, and the, the thing there is that, that origin destination matrix basically are numbers about relationships, right? So we have here, let's say we have 100 counties or 100 districts and we have origin destination data. So that means we have for each district we have the number of people moving to another district. Yeah, this is sort of the lifts, people live somewhere and commute to some other place and use a particular transportation means. We have for every transportation means we have the number from every county to every other county. Yeah, so we have sort of a matrix, right? A table layout of the origin, the origin regions and the destinations and then the numbers are numbered people. And then we have off this table, we have several, we have the numbers where people take the bike, the numbers for people take the train, where they take the bus or where they go by foot. And so that is a two, two layer, a three way array, right? With two of the dimensions are spatial, yeah, are spatial vector. They are the, the districts. And here we see the district and district of one zone lighted, lighted up here with zone 33 because we will focus a little bit on that. And if you look at the original data, the OD data, it has the origins, it's a table, it has the origins, the destinations, it has all moving people and then sort of split up by travel mode, right? So we have, we have essentially the sum of these four variables here and then we have them split up. So we have, but they're all moving people, they're just in different modes, right? So we have, it's, you could argue this is not tidy because it's all moving people. We can't sort of make this dependent on mode because mode is sort of distributed over columns. So this is one way you could organize these data. And of course the advantage is that you only have the combinations of, you know, you don't list probably this table doesn't list all the, all the zero entries because there might be a lot of zeros, but lists all the positive entries. But as you can, can guess now, we are, we're going to reorganize that. And so we do so, we're going to do some, some, some gathering. So we do this with tidy verse and here's the gather and spread, gather and spread are there they'll call pivot wider and pivot longer. What Godter does is basically that it gathers a couple of variables, sets them sort of stacks them. Yeah. So it sets the foot and the car. So we get one variable and then it repeats all the other things, the destinations and origins. So we get an OD with the mode, yeah, the bicycle and the number. So I now drop the four, this, this all variable, but I'm going to stack all these four variables into one thing. And then I need to make a couple of manipulations to make sure that I fill a matrix or an array in the right way, not going to go in all the details of that. But what I end up then with a little bit of, a little bit of magic, basically, I have to do some matching here because the, because the table is not, the table is not full. It just might constrain, it might have dropped all the zero entries. I'm going to make this into sort of a full array that has, that has two spatial dimensions here origin and destination, O and D, and it has a mode. And we see that the spatial dimensions are multi-polygons. They're identical, of course, and they are, they're geographic coordinates related to the WGS84 datum. And otherwise they don't have any information. But what you have is the values, is the mode, essentially, which is the label of the, of the, the fourth dimension, which is the, the travel mode. And then they can do all kinds of manipulations that are of interest here. This is the full table. You see that there's a lot of zeros now, but there are no NA values because of course the NA is, will be filled with this is contingency, this frequency is the NA is will be filled by, by zero. And they can plot for instance here. And what I did is do a subset and I sort of subset means the first comma means that I am going to use all attributes. The second means that I'm going to use all the values from the first dimension. And the third means I'm going to use the 33rd of the, of the second dimension. Second dimension is destination. So I'm going now to do a slice of this three dimensional cube. I'm going to slice through the, through the destination dimension and take destination 33. Right. So I get now four maps, the bicycle, foot, driver and train maps for all the people that trend that move with this to destination 33. And that's the highest, the destination with the highest sort of, it's a city center. Right. And all the people going to destination 33, where they come from, how many and what is their mode? What is their transportation mode? Right. That's in these four maps. And you can do just do that by say, give me this plot with destination 33. So I slice this two dimensional data, do a slice over one of the facial dimensions. So the other one I can map, right. You cannot sort of do the map of all the, all the combinations. That would be, that would be a lot of maps. Now we can also sum, for instance, this thing over the second dimensions of the destination and then look at which is, which district is my maximum that is this, this 33. Yeah. That's the reason why it picks 33. It's a city center. And also look at, for instance, the sums of the four transportation modes, right. Sum over the four transportation mode and get this array where the fourth dimension is gone. I can also do the origins, the origin totals by mode. You see origin totals where I, where essentially summed over the first and the third dimension and destination totals by mode. My summing over the second and third dimension. So you get a lot of summary statistics, essentially a very simple code by taking these origins or destinations as a dimension on a, on a queue. Right. So here I took these sums and over the origins in the destination and implanted them side by side using a log scale for the, for the variables. So you see some kind of reason of variability expressed here. You see the destinations. You see that the origins are much more spread and the destinations are much less spread and you can do some kind of for computations on this, on this kind of things. Yeah. So this is a, an example where you, where you have more complicated data set. And again, you can do all these manipulations by manipulating the tables using tidy, using the player type of things. But it's, it's a pretty much high likelihood that it will be very clumsy and look much more complicated to, to do it that way. So that was a more complicated example with two, two spatial dimension, two fully spatial dimensions, right? So both the origin and the destinations are two dimensional spatial are, are about polypons. So these are factor data cubes. And before I go to roster data cubes, I want to sort of wrap up that factor data cubes basically could be stored as data frames or tables. The thing is that when you do with, when you work with arrays, it is immediately clear where you have your missing values or where you have your zero frequencies to see the account data. So the table that, that, that Robin provided, it's not immediately clear that there are not duplicate records in it. And it's not immediately where the records with zero frequencies are because it's basically it just stores the positive numbers. Yeah. That is that, that's the information you have. So, so that is, that's information that's convenient to fill that basically to complete that to fill that up. Duplicate records is a, is a potentially a problem. Also if you work with the same data in tabular form, then there is no index to data, right? So the order of the records is not given. And probably you're going to do a particular kind of analysis where the order of the records doesn't matter. That means there are no indexes. If you're unfamiliar with sort of anyone who ever worked with really sort of larger data sets where you, where you have to wait for an hour or days or something like that on computations and it is about, it's about databases and finding data, then everyone will know that, that having indexes in your data is, is the thing to speed computations on. Indexes mean that you know in advance, when you pre-compute where a particular record is. And if your data are in an arbitrary order in a, in a very large table, then you need to go through it and go through it and go through it, right? This is one of the reasons why data.table, the package is so much faster because it's pre-computed and stores indexes on variables. Yeah. It does basically what a database also can do and that our data friends do not do. And if you have large, that makes a lot of sense. Erase are just there. The only reason they are there is that they are the indexing. If you know that sort of, you know, if you want to know sort of the combination of the 22nd origin and then the 54th destination and then the travel mode foot, then you know in advance exactly where to look up that value. That's the way arrays work. They're not just tables. They are clever structures that where you have pre-computed where everything is. You have that in it. You have it. It's a very fast way of getting to the data that you want. And the other thing is that you can directly operate over dimensions. You don't have to sort of go through all the records. No, you just directly operate on the dimensions. Say I want to sort of have the ratio of the disease cases over the total population. And that's just two labels or something like that, right? And also, factor data cubes arise naturally when you sample raster data cubes. Sure, raster data cubes are up next. Quite often you want to sort of look at what happens in this raster data at these set of sites or what's going on for these. How do I summarize these raster data for this set of polygons because somebody wants to know this at a polygon level because it needs to be connected to social economic data or something like that. Report it. And so factor data cubes arise naturally when you do sampling or aggregation of raster data cubes to points or to polygons or even to lines. Why would you not do it? Well, tables may be more efficient. So the data cubes that I look now at are all dense data cubes where I store all the combinations. Right? A clever way would be in the audience destination case that it would only store the nonzero values, right? And that would be a sparse array which I didn't look at. And you can do sparse arrays, but it's just more specialized. We often don't do that. Things on this quite often anyway compressed. One of the reasons not to do it is rather than be smart, use raw power. So there is a trade off between the amount of time you're going to put into developing a dedicated analysis that's very fast. So the analysis would get faster and faster, but the time you spend on getting this done gets more and more. And so the total cost of your study is the combination of the time you need to wait for the results and the time you had to develop the solution. So when I talk about raw power, then you can say, well, I know how to do this with tables. I'm going to throw it at some infrastructure that is going to do this for this big table that's going to work. So this is essentially the higher idea behind Google BigQuery GIS, which is a post-CIS and steroids. So you could throw terabytes and petabytes of data against it. At very small compute costs and get your thing done on a table or form. It's just a way to work. This is also how the company Cartow essentially does everything with vector data. They're doing everything with vector data, look at their walks and so on. They're doing raster analysis with vector representations of their data. They're just either the points or small polygons, whatever it doesn't care, just because it's easier to do everything in the vector context and it doesn't matter when it is there, it's there. So you get things done. So I wouldn't say that this is the future, but it is definitely something that might be, in certain cases, the cheaper option. Then we move on to, so any questions about vector data, cubes or my remarks about big computations and so on? Yeah, there is no, there is, you know, there's not a truth. There is different paths lead to, hopefully, two identical results. They should. Okay, then we'll move on, yes? Yeah, these great global grids are an interesting topic. So I'm not really going to. There's a person who made a package and now he's going to present the. Richard Barnes developed a package called the DG grid, this great global grids for hexagon grids, the covering the earth. Together with Dewey Dunnington, I did a lot of work the last couple of months on a package called S2, which has a different discrete global grid, the S2 grid that powers all the Google infrastructure, which is based on a cube representation and sort of quad trees on cube sites. And these are essential indexes on spheres. So when I talk about array indexes, the typical raster thing that I will now go into is like, is it two dimensional? Is the idea is that you have two dimensional space and have an image with rows and columns and that that is fast. But if you have mobile data, then the flat representation is not typically not so useful. You could do a longitude latitude grid of a certain grid size, but then your cell sizes along the pole are very small and they all touch each other. So you get these problems. So discrete global grids are efficient ways of cutting up the earth, the sphere into similarly sized parts that can be hexagonal, H3 and so on. And they can be nested or they could be semi-square like the S2. Yeah, it depends on the problem. So if you have global data, it's not a bad idea, but it really depends, I think, on how important it is that your cells are reasonably equally sized. If you do computations anyway and sort of, if you're able to compute the size of your cells and take that into account, then you're also fine. Right? Then your cells are just un-niggly. So it's not like forbidden to do longitude latitude grids. It's a stupid thing, it's just very inconvenient if you're at a pole, if your region covers the pole, it's extremely inconvenient. So then people sort of do that and it turns the earth so that the equator goes over the pole and that kind of strikes. Yeah, so that's somewhat, you know, it's a very interesting discussion and irrelevant and it's a way, there are very clever ways of indexing spheres, indexing the globe that are good and hexagons is one of them. You could do it this way and you just have a list of hexagons and put it in an SF structure and it would work with stars. It's just that stars doesn't know about the indexing very well. It just says, oh, this is just a very large collection of hexagonal polygons. Whatever I'm going to do with that. But then if it would do geometrical operations on it, it uses S2, which is again, blazingly fast, right, that has a different indexing structure. So raster data cubes. One more question here. Yeah, sorry. Can you repeat the question because of people online, people complain that they don't know about the data. Okay, I'll repeat the question for the back end, for the online. Yes. So stars and information for spatial temporal arrays. Arrays. Arrays. Yes. So, it's not necessarily spatial. It's just in the beginning. It's only non-spatial. True. Yes. It's just we want the dimensional. Yeah. Yeah. It's so, right, so repeat. The question was stars. What does it stand for spatial temporal arrays? But the examples I started with are essentially non-spatial, non-temporal arrays. They're just generic array data and that is covered as well. It is just that it doesn't have that many benefits in that case. Yeah, if I just have labeled arrays, I could work with raw arrays and do everything with apply on them. Yeah, so they then also the label, the dimension labels are being kept apply sort of propagates the dimension labels. So our has a rich sort of infrastructure for dealing and handling, subsetting, slicing, analyzing, erasing, reducing arrays, permuting them, combining them. So in the package A bind that stars also uses. Yeah, but that and on the non-spatial, on temporal things, stars adds very little. Yeah, there's little benefit to using it. It just sort of passes it on to what are already does again. So raster data. So here is a test.drd file that I read from the raster package and I read it with the command read underscore stars. So read underscore stars twice to read raster data using GDAL. And this is a vector format that actually is dedicated for our raster data sets that GDAL can read. And that we get here. So we get the sort of the test.drd variable. You see that I get two dimensions here, x and y. So I now have two spatial dimensions, yeah, two dimensions, x and y that kind of together sort of describe space. So this is a very extremely different situation. I have x and y and I need both of them to sort of find a certain location. So in vector data tubes, I had just a sequence of either points or polygons or lines, whatever, that each of them described two or three dimensional space. Here I need two dimensions to essentially to get that. I could select a value sort of for x and take x with value three and then I would get a constant x meaning I would get a column of my raster out of this. I could select a column or I could get a column subsection of that raster by selecting x values. So we see here that several things. We see that x ranges from 1 to 80. So it's an 80 by 150 array, but that we also have offset and delta values. And offset means that is where the raster starts on the x side. So the x edge of the raster is the offset and the y edge is the y offset. And since y values, y delta, so delta specifies the step size, so the grid cell size in x and the y dimension, you see negative y values, negative cell sizes in the y dimension for the reason that the y coordinate increases when the row index goes down. So they have opposite directions. So that means a negative cell size. Why goes up going north? The typical images have cells that the image, when you plot the image, it goes from top to bottom. This is the way image data is almost always stored. You don't have to do that. It can also be the other way around, but almost everything that is stored in file is stored that way because PNG, JPEG, everything does that. Images start bottom to top. So that is the delta. And we have a reference system, which you can see here is this 4 string, which are now deprecated. And the other things are missing. We have a special label that says, is this one of the x and y? Yes, this is the x and this is the y. The reason is that not always x and y are called x and y. They might be called latitude and longitude or lat and long and so on. And somehow they meet, they're mapped to what is going to be plotted as the x and what's going to be plotted as the y. This sounds trivial and it's not always trivial. There are a lot of cases where if you have a polar projection, it's not clear what is worth it. And then if you plot it, we use this theory, this scale, and we use equal breaks and we get the usual plot of this new data set. So we see here a couple of things. Raster files are data cubes. The attribute has the values for all combinations. So also the white cells here are filled with NA values. You see here that I have 6000 NA values. So this is all sort of the whole roster is filled and I just plot, I just show the colors for the non-NA values. Their vector data cubes have a single spatial dimensions. Raster data cubes have two spatial dimensions called with x and y labels. Or they might have three, if it's a three-dimensional data cube with foxholes or some kind of meteorological data and so on. This case is an example of a regular grid. That means that offset and delta are filled and that indicates that we have constant cell sizes in each of the directions. Here there are also square cells, but as you see they don't have to be square cells. This one might be 40, this one might be minus 50. It's fine, but they're constant over the space. That means that I can compute my coordinates from the index by this equation. If I have offset and delta and I want to know the coordinate of index i, then I can compute it this way when my index is one-based and then for the value one, I get the cell edge and for the value 1.5, I get the cell center. We have the cell space, which is continuous, and the index, which is also then thought of as continuous. This took a while before I understood. Then why is the delta for y is negative? That means that the coordinate values decrease when my index increases, when I'm going down. Also, as I mentioned, you have a reference system. That means that I can link this data to other spatial data and put it in some place in the earth. What is plotted are the non-missing value cells. What it can do with this data is convert them into an SF object. That means make them vector data. What then happens is that, and I say, don't remove missing value, which would by default happen just to show you. I get polygons. I get here then 9200, which will be the product of 80 and 115. Polygons, and if I plot these, then I get a plot like this. You see the little polygons plotted, where I also plotted the polygon borders in gray. You see all the grid cells plotted as small polygons with a gray border, and the values filled whenever there are non-many values. You get the same map in a vector representation. You can equally do that. This is what I meant by what Carto does and what BigQuery GIS would do. Look at the whole thing as a bunch of polygon data. It works also. Of course, you wouldn't do this with massive grids. This works nicely because it's relatively small. If you're going to do this with very large grids, then you're going to set your computer on fire. Not literally, but it will require so much memory that your process is breakdown. These things actually become really interesting when I use multiple layers, raster data cube, so-called raster stacks. Yes, your question? Just a brief question to something. You said x i, where you indicated x i with the offset grid. So that's the net edge. There is a value offset and delta that you find in this table, offset delta. You filled it in in this equation, and you give it an index. Sure. I know about that. That's the net edge of a box, say, area. And why? Is it the top reference of the box? The question, repeat the question. For x, I can see how this works. It's the left side. For y, is this the top or the bottom? And here is the magic. It depends on delta. Delta is either positive. If delta is positive, then an increase in index means you get an increase in y. You go up. And the offset will be the bottom edge. If the delta is negative, taking an index of one, the first thing, cancels this element, means that the offset is the top edge. So the meaning of offset depends on the sign of delta. That is the trick, how to think of offset and delta. Thank you. Good question. Right, multiple layers. So here's an example from the RGDAL package, which has a familiar JPEG image. You can also read JPEG images. Reading JPEG images is not a problem. You get three bands, RGB bands, because it's a JPEG. And you see that here are dimensions. And you see that these are regular discretized. I again get an offset, which is 175, and there's a delta of minus one. So again, it's going down, indexing. But the reference system is missing. So it's regular, but this is not Earth data. This is just an image. And that is also why the cell size is one. But JPEG can read such things, nevertheless. But it's just unruns you in spatially referenced. You get the intensities in values from 0 to 255, which is typical for JPEGs for the three bands, R, G, and B. And then you can use the function SD and RGB that computes RGB values. So R has values for colors. You can use the RGB function to compute a, ah, it needs to be maxColor is 255. It gives you labels. This is text data. This is character data. It gives a hex representation of a color that can also include an alpha value. And this basically throws these three bands together, or possibly an alpha band as well, and converts that into color values. So that reduces the third dimension of the band. And I get these color values out. And when I plot that, you see the regular RGB thing. I don't know why it's not plotted. It's the thing you've seen before. Let me try. Here it is. So I need to. I'm not so good with. I'm not so good with our studio, to be honest. This is very embarrassing. Oh, yeah, this is it. Yeah, thank you. Yeah. Thanks. And then you get the unsurprising, you get the color map of the R logo. And that puts RGB on top on computer colors. Right. So this is something you can do to compute colors from values. Right. Multi-band data is more common with your handle satellite imagery. Here's an example of an excerpt of a small Lanzon 7 image used to cover part of Brazil. You see it has six bands. It has all the 30 meter bands. And they're there in UTM zone 25. Yeah, so Lanzon data is in UTM. You also see here that the point label is set to false and no longer to N a. That means that the data set that we have been reading here that this stiff file has metadata. That tells you that pixels are actually the pixel values are associated with special areas and you're not associated to point values. Yeah, there's a difference between one and the other. You could have other data sets like, like the M data digital elevation model data where your pixel values refer to an elevation at the cell center. That is a different that has a different semantics. If you then would sort of zoom in, you would like to interpolate values rather than draw them as square boxes. Yes, here you hear interpolation is sort of harder to imagine. Of course, this doesn't mean that the cell values are a constant value throughout the cell. They're just associated to the square and they're of course they're confolutions because this isn't this is an image sensor that that sort of measures averages over small areas that are not even square. Right. And if I plot that, then you get you get a plot like this, and it essentially histogram stretches over these seven bands, right. So it takes all the values of these of these array. So this is a 350 by 350 by six array takes all the values, and then takes the quantiles and sort of assigns course that means that for these six maps, the sort of the 10 different grade scales are kind of equally distributed. Yes, so it's similar to what you do with high dynamic range photography that you kind of distribute colors equally not get very large, light or dark areas that you see something if you would do it in in standard constant or equal interval range the way that roster does it and most of the times all all the maps will have the same value will be all dark or light you don't see anything. Yeah, so you can override that here I did the same thing but then I didn't do this with a joint Z limit so here all the maps sort of refer to a single color scale but I can also do the histogram stretching for each of the bands. Yeah, so I do that here so here you see more color stretching and of course no key anymore because that each subplot now has a has its own key and has its own color breaks. So you can do this to sort of blow up variability. And it also maximizes spacing of bottom region so you did you see something right instead of putting legends everywhere. Right, so you can then also do RGB composites here this is the true color composite this is some kind of false color composite by by computing by computing RGB this is by plot and using argument RGB I could also do as the RGB with the band selection, and then compute the RGB was different ways of doing this. Transforming and warping rusters. There is the thing that if I want to sort of combine Sentinel to and then some eight and it might be the case that these two don't have the same code and reference system is definitely the case if you put something on top of modus data which is this sign is vital grids, then you need to work with that. You need to do something and this to illustrate that I created a little data set here. This is this bounding box. Yeah, this is pretty much sort of a large section of Europe. And then I said, converted into a stars object, which has essentially an empty at has zero values. And then is a 30 by 20 grid with with one degree grid cells associated to WTS for covering that part of Europe. And then I did this not going to talk about what that means. Maybe in break and loaded some our natural Earth data, which is which is global Earth data with some, some, some, some parameters, some parameters and downloaded the Europe, the Europe section, which includes Russia, and then computed population densities by basically taking the population estimate and dividing that over the area of each of the of each of the polygons this is polygon data. And then I convert the units to one over squared kilometers yet to not get very crazy values, very crazy low values because the this is one per square meter. So I put all these data, and I'm next going to rasterize that that means that I'm going to take the population density polygon data and rasterize this over this over this object, it's that they just made which is the one by one degree grid. And it just sort of takes the does the population density values for these grid cells, right. And I do that. So it's very course map of a long acute latitude. You see that, although the grid cells are square in terms of degrees. Yeah, you think square grid cells one degree north and one degree east. Of course, one degree north and one degree east doesn't everywhere on the world mean the same distance. Yeah, so you see now that the plot methods, given geographical coordinates taken aspect ratio so that in the map roughly one meter north equals one meter east, which is not the case with degrees. Right so it doesn't aspect ratio correction on as what is essentially an equidistant equidistant rectilinear projection. So it's not degrees let it as if it was x and y, but an aspect ratio correction to make sure that countries have the kind of the shape that you're familiar to and they're not like like Germany is this compressed term this kind of round Germany which you see when people ignore it. So, so this is something you do and if you do that then you also see that the roster is no longer has no longer square grid cells. Yeah, but has grid cells that in kilometers in the center of the map, roughly, roughly has one degree in one degree and see that one degree latitude is like point eight degrees long. That should be the other way around. Okay, so what I can do what I wanted to make the point now is if I if I'm going to transform this to another coordinate references for instance the ETS 89 Lumber equal area projection which is the recommended European sort of plotting which is at least equal area plots of Europe, then I can do this for this grid for these great values and I get something that is that now looks the same is the same data and is still a 30 by 20 grid, but has a different reference system and has a big bunch of coordinate values and it's now it says I'm a curve linear grid. Right so what I essentially did here was nothing but just change the coordinates of every grid cell. Now here I knew the coordinates because I could compute them from my regular discretization. Now I'm moving them to a new system which is either where I have to compute for every grid cell, I have to compute the coordinates in the new coordinate system, and they look like this. So it's still a regular 20 by 30 grid, but every grid cell now has its coordinates in this in this coordinate reference system which is EPSG 3035 which is the ETS 89 LAAEAEAEAEA projection that actually every European citizen should be aware of that this is the way to plot your European to use for your European maps. So it's just a little bit of a digital familiar. All European publications will be will use this. I think. Right so it's still a grid. It's just not. It's no longer regular. It's no longer regular. It's so called curvilinear grid. And that is that can be convenient or inconvenient convenient thing is that all the grid value are still there. I just changed the geometry. Yeah, I just the geometrical transformation. And of course, at some stage you want to work then with our systems that need regular grids and so on. Otherwise you can then do is warp this grid warping means that I have some kind of up shaped thing and I'm going to burn it into a new regular grid. Yeah, or you can work from anything to anything. But this is what's happening here. I take this curvilinear grid. I'm going to make a target grid. This just the sort of the default grid sizes that takes the bounding box of my population data on this reference system. And then happens is that it's going to redo the values. Yeah, so here you see then the plot of that new grid that you see here the sites, you see the ruster cells and much much much smaller. And of course in the values you still see the lucky thing of of of what I had here. It's just rusterized to a new roster which is a regular roster in ETV RS 89 L a EA. Yes. I don't know how to open the chat. Yeah, then I get back to my thing. Yes, you can load the news data set. That's a good question. So, if you don't have this data set you may have to install the roster package. Yeah, because that is where it's going to read this from this data from if you install the roster package you will not have no problems to do this. Okay, good question. So this is what warping does. Yeah, so it works one sort of roster in one corner reference system to our corner preferences. And it means it's going to recompute values. Yeah, so the totals of your roster will no longer be the totals. So for density, they still makes some sense for population density. If I had done this for population numbers, it would have been entirely senseless. So it should be a population here over roster cells, and then re rasterizing that will sort of destroy the whole idea, because they are kind of extensive variables. So, in that case, you would use as the SFST inter interpolate area way that use area weighted interpolation to go one from one, from one rasterization to another one and do that on a polygon basis. Yeah, the way Carter would do it. So now it, yeah, sorry, any questions to related to the multi layer rosters, curvilinear grids grid warping RGB. Then we will go to the more exciting things, which is roster data cubes, time and multiple attributes. So we haven't seen any time variability here. Now we're going to load a data set that is small data set it's available in the package stars. That's a net CDF data set right net CDF is a format that has been around. It was, I think it was one of the first sort of libraries that I discovered online. This is in the, in the, you know, decades before HTTP was invented is in the time that we used ftp to, to sort of discover data sets right you would ftp to some first of all you would go to, to groups to use net groups and see what's there right to read emails and sort of mailing lists distribution system and see where you could find things and there would be no outsides or something like that. And they would have net CDF files and a library called net CDF that way you could read and write these things right this thing still exists and has sort of developed over the years and is still dominantly used by modeling communities right so it is. They use it for spatial temple they use the grid of data. And, and for data where you have multi dimensional grid typically time varying things like weather data, like Earth observation data a lot of things are nets. So here, this thing is being read and you can see that it reads to sub data sets or to sort of attributes in this data set, and it is a three dimensional grid layout, where we have an X and Y. And you see that the deltas are one eighth of a of a of a degree. But you don't see a reference system. Yeah, so that's unclear this this is degrees, let it be long queue but it doesn't say so. And then we see that there are two value precipitation and something with temperature. And you see also that that it reads units. Yeah, so these data sets come with actual measurement units, which is which is very cool. And they are being read and they're being propagated using the units package, which is another story, which, you know, which is extremely, which is an extremely good idea. But occasionally causes some headaches for people because you sometimes you want to, you need to get rid of the units before you can do a next step. Yeah, if you want to put things with Gigi plot that not always works directly. Here we plot this data, and you see what it plots is a sequence of maps that loops over the last dimension. Yeah, so these are 81 by 33 roster data sets, and they have 12 of them and you'll apply them in the labels are the dates that are associated with each of the sets. And these are apparently monthly something monthly averages or monthly totals of precipitation, probably total precipitation is an average temperature. Otherwise, it would make sense. But we would have to look in the metadata of this data set to get more information about who created it for what purpose what it is exactly I think this is climate down scale climate data. Another example here is where we have data that are distributed over that over several net CDF files. And this is typically by monitoring programs where a satellite flies over and measures something and somebody realizes that data and sort of sort of converts it into some kind of product like sea surface temperature and uploads a new file every day. Right. Instead of sort of rewriting and very large pilots and larger and larger, the file is added every day and you see here that the that the date is part of the file name so that makes it easy to to find these things. Yes. Oh, thank you. Oh, that's good. Yeah. Thanks a lot. It should work. Yeah, it's half an hour to go. Other questions or warnings. I have a set of net CDF files here and it's going to read them in a single bunch and it's it the read stars here is going to sort out that all these are also data sets that seem to be arranged along a time axis it's going to go into each of the file read the time stamp for it and then sort of collate into a regular time dimension where you see that I have nine time slices and that this is the first value so this is an offset this is a post CT structure so you see the post CT value of the of the first day which is two hours Central European standard time so it's midnight UTC just my system wants to plot it in in my local time zone, which you can change and has a delta of one days. Yeah, so it sorts out that this is regular time this is a regular time axis. And it doesn't know whether this is associated to time point values or time intervals that is not in the metadata. Although net CDF actually can have that as metadata it can if you do temporal aggregation, you can say my temporal there does time bounce and sort of it's really an average or a maximum over a certain time period and then it is weaker on the on the spatial thing because it doesn't have a current reference system for whatever reason that is probably encoded in a way that is not interoperable or that she doesn't doesn't recognize. That doesn't, she doesn't recognize this as degrees in some kind of associated to some kind of date. What you also see is there's a fourth dimension which is called Z level. Well, this is C surface temperature. This is a surface right so it, they're not much interesting going on the elevation zero, it is C surface. So this is zero meter. So it has essentially a unit for the elevation, but it's a single layer elevation. So we know it has metadata. This is on a certain elevation that I know zero meters with respect to C level, but it doesn't have any variability here because it is surfaced. It causes a bit of trouble sometimes you want to plot it. It's going to loop over the third dimension and so show, see, show the sea level which is not terribly interesting. More interesting is the variability of time. So I can drop this. I can use a drop to drop singular dimensions dimensions that have only one value. It's a, it's a 1440 by 720 by one by nine array. But if I dropped it one dimension that has only one value, it remains the same information. I just lose the information that this is at zero level. So I have to, so it's there it's in the data. It's a four dimensional array, but with so to speak and obsolete third dimension. I drop drops that one and retains the other three that have variability. And then I could plot it. And here is an example where I used Gigi plot and GM stars with a couple of other things to to plot these data and the various thing and so on. So basically you, you give it, you give it Z and then see one means the first attribute. If I don't do that, it takes the first attribute anyway I have a certain alpha value, and I'm going to down sample this yeah go to down some of that means that, although I have quarter degree, quarter degree grid, I'm going to down sample this to a two and a half degree grid, because my plot is anyway so small that every single pixel is probably already five degrees or something like that. So, I don't have to do that. But I'm doing that because otherwise it takes the eternity before this plot is finished. Otherwise it's going to do all the pixels. Yeah, and then it takes like five minutes and every time I compile my slides to wait to five minutes. I don't know if it's five minutes but it's it's a factor hundreds that I'm reducing here in the number of pixel data show. And in the end, the thing looks the same. Pretty much the down sampling is something that we haven't seen before, but it also happens here it happens in every plot. So, I have a lot of different commands, but these are base plot commands and there I can, you know, when I do a base plot can check what it was my device how many pixels do I have how many pixels is area have here I mean ggplot and you can no longer do this because you are sort of lost in the ggplot abstractions of everything. And at the moment that ggplot, you know, ggplot distributes tasks at the moment that you're in a geomunders core something you cannot sort of ask ggplot how many pixels does my device actually have. You cannot do that. Yeah, interestingly, I had, I had discussions with some of the real sort of ggplot gurus. Thomas and Claus Wilka Thomas better than so on and you can't you can't do this. And it's very, which is very regrettable that you can't do this. Yeah, so sometimes ggplot will take eternity to finish. Yes, and then remember that you can do this down sampling by hand, regrettably, which you could do automatically if you knew if if if you know that. So this is developer talk, but also warnings that some things, you know, are difficult to do in our other operations on data cubes are subsetting we already saw if I take the first sort of a single index. Then what happens is I just do a listing. Starsh objects are lists with a class in a dimension table in an attribute. There are lists of arrays. So if I do a selection of the first element, I just get the first list of the second list element which is here, the anomaly. I just drop all the other arrays and keeps the dimension table because that doesn't change. Other things is that I can select other elements by saying comma comma means that I select all the first attributes. I don't specify anything I select everything. So it selects attributes is the first square bracket selector, then it selects first dimension second and so on. So first and second dimension x and y. So now I selected a subset that has the first 10 rows in the first 12 columns. And then it's reversed the first 10 values of the x variable meaning the first 10 columns, and then the first 12 of the second dimension was the first 12 rows columns rows. This is how the thing is organized. So there is this, there is this fortune about illogical, which is, which is one from Brian Ripley in 2004. The medical conventions are just that conventions they differ by field of mathematics. Don't ask us why matrix rows are numbered down but graphs are numbered up the y axis, nor why x comes before why but row before column. Right. This is what I just, what I just mixed up. Yeah, and you always a row column, and you say XY, but it's calling wise row. Right. So if you talk to me this, this is just conventions. So I see level is still in there and retainer. So, and you can also select and see level of course they only have one there, and I can also select time. Yeah, that means that I need four commas right for the attributes x y and z level, and then I get a smaller time slice right I select from time time item five three from three to five. And now you see why the from is there, because I selected from three to five, but I still have the same offset. Yes, and if I would say okay, those are three time instances, I do from one to three, then I need to re compute the offset. Yeah, because the offset shifts. And that's why this is why this from field and this is sometimes confusing there isn't st underscore normalized that that does that. Yeah, that will sort of change this into from one to three, and then modify the offset value with the appropriate number of deltas. So, so it can also be used to shortcut intersections as like here, if I do a global grid, and then I plot sort of this, the section that intersects with that intersects with the Europe, European polygons. Let's see if we can do that. I think I didn't do that because it looks very ugly. And it looks very ugly, because, because any is is European data set that includes entire Russia and that crosses the deadline. So you end up with a grid anyway that is global because the deadline is involved and all the kind of I didn't actually, I didn't actually do this. But the idea is there right so you have you do an intersection you select by, by sort of the area you basically crop the area by a certain geographic region that can be a simple feature geometry. So you do is cropping or making sort of a sub sub roster with smaller area, or you crop by a polygon set, and you then basically and burning and a values in all the cells that are outside this area. So you have the verbs that tiny study verse type verbs that work with stars object. There is a slice that selects dimension values index. There is filter that selects dimension value using the dimension value so you filter you use coordinates for slice you use index you use there was pools that pulls out an attribute select that selects attributes and mutated computes new attributes based on existing ones. Yeah, so you could mutate and then say compute sst times five or something like that or some conversion into Calvin or some conversion into far and height or something like that do that on an area based basis right for all of the cells. Then there are also kind of these kind of map algebra type things we can load in the set from from the from the launch out imagery, and again, for instance, multiply by two, and it gets something that go to 255 because to 510, and it can also do things like mask things out so so for every for every cell that is smaller than 125 give you a value in a or something like that, and then do is basically do these kind of array array logic right so you basically retain all the cells that have values large under 25. So a number of sort of things that that that you could do that you think you could do with these types of things like algebra like, like expressions you could actually do with these stars. So I just saw the sdssf, which basically takes rasters and then draws either points on all the roster cell centers, or draws little square polygons around all the rosters, depending what you want. That is one way. There is the raster eyes that we saw that goes from vector to roster and raster eyes is then that essentially calls the GDAL utility rust. It's not. It's part of the GDAL library, the raster eyes routines there. And then there's the issue of proxy objects. Yeah, so and so there's a couple of ways of going from raster to vector and vice versa. There's also aggregate cells over polygons and computes whatever function to some or the mean over that polygon of all the values. There's even a very fancy one that does that using the package exact extractor which is a terribly fast and good package for for computing aggregates over polygons that basically looks that cuts all the little square raster cells into parts and sort of use the part of the raster cell that falls into the polygon part of false hours. So you get this really, really area weighted aggregations for some some means. And that does is an incredibly fast way with a very clever algorithm. Exact extractor. Exact extractor by Daniel best one of the SF contributors and Geos developers. Then there's the thing how do I how do I work with large raster data. Yeah, so I have my hard drive is full with Sentinel files and so I'm going to go through that package raster does that and in in a particular way that is not not too transparent. We do that in stars with do that in somewhat different way in sense that we have this idea of proxy objects. And I'm doing that here with with Sentinel to say that you can use a granular or tile or seen or whatever you want to call it that is available in the star state of package. So I'm going to create here this file name which is a very long file name which has a prefix, essentially to tell GDAL this is a Sentinel to L one C data set and tell you where it where it is and then give sort of the long thing of it need the thing that it needs to read inside that zip file to get a metadata, and then it reads it. Yeah, so this is all, you know, this is all something you need to do once and then for every time it works the same thing. So it's going to read it here. And what you see is that when you get back is the stars proxy object with one attribute in this file. And this has the four 10 meter bands. You see here that bands are named before the five p six and the eight I think is the one that I have 11,000 by 11,000 pixels, each of 10 meter size, and they're in a particular UTM zone. So this is the four 10 meter bands in a Sentinel to tile. And it's not reading anything. So this takes like less than a second because it doesn't read anything. It just says the data are here. I know where the data are. And this is the metadata dimensions. Yeah, so it tells me sort of the dimension of this area without reading it because reading it would would sort of take a lot of take most of your memory. We can control whether this is being read as a star's proxy by setting proxy is true but at current state of affairs. It's being read as a proxy object by default because it's so large. So it's first going to read the dimension of the thing. And then besides, okay, there are more than 10 to the power eight cells. I'm not going to read this into memory because most likely you don't want that. I'm going to take the have this as a proxy object. And then it tries to do lazy evaluation that means that it doesn't do anything up to the moment that you wanted to do something. Right. So at some stage you want to do so you so you can compute for instance here. And the VIs by applying a function called any vi that called that computes your any vi values over the band dimension to the accent to look that over the X and Y dimensions right so apply it to the bands. This is this is the expression so band four is near and then one is read. So if you do this, it's again it's done within, you know, within a microsecond because still it doesn't do anything because you didn't ask for pixels. It just says okay I still have this property object. It just now has a call list I need to do something. So if you want to have something and going to do something. Now if you plot this thing. It also does that in like something like seconds or so. Because because what it does is basically it computes as many pixels as you have computers many cells as if you have pixels on your screen or on your device. Right so you can do this on a 10 meter resolution computer and the vi's and it would take like 10 minutes to read the whole thing and then compute it and probably write it back because it's too large to hold the memory and so on. So it would take 10 minutes but if you only want to plot this then you basically look at the plot of what will it be 600 by 600 pixels or something like that. And I'm now going to loop over the 600 by 600 pixels read them compute any vi's and plot them. Right so this is a bit of a trick so so this plot that any vi actually looks at this object and say okay there isn't this object but it still has homework to do, but I need to plot. So let's look at what the homework does. Oh it's an sd apply or it doesn't do anything special it only works on the bands so I can do it for just as many pixels as I need is no filtering going on or if it was filtering I of course couldn't do this trick. Yeah so it looks and says okay I have to do something and do this at the screen resolution and then plot it and that speeds up things from, as I said from something like 10 minutes when you would do it for the entire day to something that this within seconds done. And the trick it does of course is that it uses it in which permits is usually over views this is rush to tile comes with lower resolution overviews of all of the bands. And if you asking sort of a down sampling of whatever effect of 600 and it can use your work on one of these overviews and only read a very small amount of the actual data set. This is a pattern that everyone who works Google Earth engine is familiar with. Basically, you put it, but that does is that it just computes the screen pixels. And then you zoom in, it's going to read compute the screen pixels. So although you work on massive imagery, you just do computation that are sort of proportional to to the plot size you want. And finally, a more complex, a more real example that that I sort of cooked up in the last few days because I discovered that Mario's apples material for tomorrow was was online. So I know I'm always curious like what's going on there. So my years will do a sort of a serious sort of quasi serious scaled data cubing. And it's not going to work with petabytes, because you don't have petabytes on your on your laptop, but it's going to work with a sets the size that you can still download it and work locally to get the idea. And Mario's compiled a number of lots of eight scenes over for a particular area over a particular time range and creates cubes from that. And I'm going to give a very poor man's poor man's version of that. Yeah, so what I do, for instance, is I look at only at a single scene so I can I don't have to, I don't have to worry about all different places and times and so on that typically satellite scenes they're not they're not in the cube right so it's an orbiting satellite so it's a trajectory of images you can. The stars hasn't as the mosaic so you can mosaic things together but you then need to first sort out which things to mosaic. And that is exactly what G dot cubes does what I sort of did not do here. But then I downloaded Mario's data set, which is, which is here you can you saw all these all these tiles here which are in directories and in in one of these directories you see sort of you see a number of tiff files. I'm going to. And we're going to load these seven bands of these of these duties. And do this in this particular scene. So, I'm going to match all the file names that that match to this patterns or band one to seven from essentially from this, the directory that starts with this expression. And then you read all this so that gives me like 120 files. I'm going to read all of them in memory so not as proxies and memory. And then I'm going to warp them I had to do that because there was tiny differences in there couldn't sort of put them as a pile on top of each other. And then I'm going to do a slight differences in in the number of rows and columns so I kind of want them anyway to to the geometry of the first one right hopefully not too much lost. And then I'm going to glue them together which can be done by the sea basically I do a sea with with the 120 objects, combine them and and then I get essentially that that set. And then I got the time series actually from the scenes somehow right because the time. Yes, you can see here, the time is encoded in the scene name. So this is going to happen here I'm extracting this, this little section of the file name, and I'm going to read that with this format string, and I get my times. Sort of obscure way of doing that. I wouldn't recommend it to anyone, but it's possible that then I'm going to take the band name from the, from the following that takes this little string. Right that tells me the band name. So I get band 127 and so on. And I'm going to merge this whole thing. And then I see that I have a seven 60 by 70 70 roster with 182 layers. And then I have to assign band and time because this is these are the two layers are seven bands for like something like like 40 different time instance. So I need to then read and mention this thing. And what I sort of do and now comes the ugly part. So, so this is why, you know why you need proper packages that do this in a sort of in a readable way. I'm going to sort of just work on this raw dimension table, remove this attribute and then create two new dimensions. Actually, one with the unique dates and one with the unique bands, and I'm going to read and mention this thing. And then you see the dimension table looks like this. So I have now 26 dates. Yeah, which has all date which apparently have no offset and delta so they're not regular. They're not found to be regularly distributed time. And you have seven different bands. And the question is whether everything is correct now and I'm not entirely sure whether everything is correct. I just did this. No, everything is correct. This might look okay so here I do a slice over date so I take one particular date. So the first date which will be the 6th of August, and then you get a plot of the seven different bands for that date so this time slice. Now this was the essential him thanks by that they did all values more than zero were missing. So I had to assign them an a this this you can do this for a four dimensional area like this. And here I did a plot of a time series of the first band. So this is band one and a time series and I'm entirely unsure whether this makes sense. But anyway, you can you can get an idea time slices, then slices, and look at these data as a data cube and sort of get through them. And you can even compute NDV eyes by doing the same trick with them with then five and four which apparently correspond to and near infrared and red and apply this to the first three dimensions x, y and time. So we compute any vi over the we use any vi function to reuse the band dimension and we get back three dimensional data cube. It is not just time series where where band was reduced to end the vi. And we get these in the vi time series, not sure whether they really make sense. But this is what what came out. So it's still kind of work in progress. But just to show how you can manipulate these sort of. This is all in memory. But in principle, you could do all these kinds of things out of memory, but and sort of still have plots exploratory plots and do simple functions on dimensions. In a relatively responsive fashion. So without doing the full computation on the on the entire. I hide dimension data cube, right, you see very small surplus. So the number of pixels that are being computed is here is also very limited to just kind of the same number of pixels that you see here maybe 100 by 100, rather than the rather than sort of the 700 by 700. So that is a kind of a factor 50 speed up that you get in the plotting, plotting action. So these were some kind of, you know, some demonstrations of how you can juggle around with with factor data tubes with raster data cubes, where it really becomes exciting if you have more than four dimensions, more than three dimensions. So more than three dimensions is also the, the use case where, at least as I understand it, the roster package and its successor. And kind of where you get stuck, because it's, it has a roster stack, as it's, as it's multi dimensional thing and if you have a stack of bands and times, and you need to do the whole administration like where are my times where are my bands, and you can you can do a function through that. It's not two dimensions, right. And the data structure you can handle is really three dimensional doesn't sort of work at more dimensions than that. So that was one of the motivations to do this work. So are there any questions from your side, which I can take and then repeat so nobody in the field will hear it. Yeah, for example, very efficient overlay for the hard to do it. Well, as I said, if you do, if you have a complex data cube object, like, like here, like in this last example. So, when I do this, when I do this plot, so if I, if I use stars proxy object so here. At this moment everything will be computed because the thing is in memory, right the whole thing fits in memory it is in memory. But if it starts proxy object, then it has this lazy affiliation so the object is on this might be massive, but it looks at the order of operations and says okay you want to plot plots how large my plotting device how many sub maps and I'm going to plot how many pixels do I need for a single subplot. Oh, that's 100 by 100. So I'm going to sub sample and just read at a lower resolution. Yeah, and that goes to GDAL and GDAL is very fast and doing lower resolution operations and you don't have to tell it. It does this because it knows it's a plot I only need so many pixels. And does the star the same thing. I also don't know there's very little information about what they're actually doesn't really doesn't do. I don't know what overlay means. Extract. Well, you know, there's a lot of things here are also done in GDAL which is C++ and that's what they're also does does calls these C++ calls C to retains. I hope it does, because we understand what they do. But it also does a lot of things in ways that are sort of hard to track down the exact details. Well, this is not I mean it's just unfriendly to say I'm very sort of, I'm a big admirer of Robert Highlands works. We looked at RGB composites of plots, and I looked at what the map view code did. To do this, and I looked at some other package that also did that, and they all would sort of look at RGB values and then take this, this two and the 98 quantile and throw everything outside of that away, basically set it to 289 quantile and then use that range to compute RGB values to do a stretching. And that happens to come out of raster code, which is not documented way you can't control the quantiles, but it does that. And it doesn't show up in the help page of the plot RGB commands. So plot RGB does that and if you want to reproduce with plot RGB does you can look in the code. You see that it does that. It doesn't explain that it does that and doesn't explain why it does that. So these are the kind of the little things that that kind of start to live their own life and then people copy that because they think what RGB is the thing that is true, right, because that was their first and everyone uses that. But I don't know why it picks the second and 98 quantile and not the five and 95 or the one in 99 or the whole range. Yeah, there are always like funny values at the tails in a lot of satellite imagery, funny high value so if you take the data range as a default which might also make sense then you get completely black images. And so there is there is stuff that you ran into which is kind of came up for that probably Robert got his inspiration from somebody else, but we don't know and it's not documented. At least they couldn't find it. So that's the essence. I don't, I don't know what all the functionality of the roster packages. That is those are two common use cases that is what it does. Yes. So if you look at sort of the, if you look at the amount of functions that the roster package provides it's, it's massive. A lot of multi core things. I think the, I think that as the apply has a couple of options, because of lies, always like in independent job because you work on dimensions. The chunks are always independent. So, as the apply has had options where you can pass a cluster and tell whether to use future or to use water parallel apply or something like that. So, I guess parallelization options. So, you know, you can use the roster object calling vocal on the roster object and converting back so that was like a three liner or so that you could throw in sd apply because the focal is essentially over x and y dimension, you leave everything else intact. So, for each of the time for each of the bands you do the focal on on the XY slice. Yeah. So if you put, I think I'm, I think I'm offline, because I did some VPN stuff that they shouldn't have done. So if I look at stars issues and then search for focal, focal filled with stars. Here's an example. So you do an sd apply with with the focal function. Yeah, that doesn't. So, the focal and roster. So does the roster conversion the fly and cost. Yeah, we could port this code is this is also kind of kind of the same thing the question is whether how this will work with with massive with massive grades. But I don't know how many people want to do focal things and massive grids. I think you can do star proxy things with roster as well. I have to look at it, because it's essentially rest of the same thing right so says here's the file name, and then I'm going to work with that and sort of return when I'm done. So yeah that is that is relatively they work to sort of stars provides the from roster to roster conversions. So you can do call call things there. It's convenient. In any case you don't have to sort of work out your band and time things. Yeah, and the other thing is the conversion integration with SF, which is pretty natural because it's kind of that sense it starts extends SF for multi dimensional feature or feature time series, or multi dimensional feature data. And then another questions. So, yeah, so, as I said, there is this arm art down so I hope that some people are some people who actually managed to run the whole thing. So in the meantime, yeah, good. Okay. So in the meantime while I was talking it's sort of the run all button on our studio after installing a few bigger packages. So when you don't manage then, you know, catch up with me in some breaks sometime. I'll be here until Thursday morning. I'm happy to discuss things are happy to also explore all the things you cannot do with with stars.
 Anyway, in this GitHub page, you will find the Rmarkdown file.  And the Rmarkdown file, you can download, and then you could go to download it somewhere.  And then you could go to the place where you downloaded, download the whole repository,  like I have it here.  Oops, a lot of things there.  But there's also the R&D file.  And you can open that with RStudio, and then you could theoretically go through the code  sections and basically run the whole thing.  So you can start there, and then whenever there's a code section, you can run this,  and then go to the next code session and so on.  And run all these things.  And hopefully this command might also work.  This command actually downloads a somewhat larger data package with some real remote  sensing image in it.  So that's a package that has a size of about one gigabyte that you can throw away afterwards.  But that I use because Chrome packages can only have a certain size, right?  So you're not supposed to send your one gigabyte packages to Chrome and have it distributed  there.  So I set my own repository where you can install things.  I think it also works with HTTPS if you prefer.  Usually you can do it just now.  Although Tom says not do it, but it typically it works, but just goes a little bit slower  if multiple people do it at the same time.  So that is where the bigger examples come from, but I will go through a lot of smaller  examples primarily.  First, we'll say something.  So what are data cubes?  So we looked at a little bit that this morning and there are essentially rectangular two-way  tables or multi-way tables where there's not a structural difference in the data.  So we think of data frames also as tables.  The data frames are essentially collections of records, right?  Like a data frame table.  The way they are stored is that they're collections of columns, right?  Because columns are vectors.  Vectors can be of different types.  So you can have a vector of numeric, you can have a vector of date, you can have a vector  of logical, you can have a vector of factor, right?  The categorical variable and essentially have a collection of these.  We think of them as collection of records.  The way they're stored is that they're a column store format.  But in essence, it's not a two-way layout table with the same thing along rows and columns.  It's if there's heterogeneity along the columns and there's homogeneity along the records.  And that is basically why tables or the data frames are thought of as one-dimensional structures.  So one-dimensional tables are also one-dimensional data.  Of course, the sort of the analogy of a cube fails there.  But if you think arrays is being generic, having one or more dimensions, it is a, the  table is a one-dimensional cube.  So typical data cubes is that they have two or more dimensions.  Then it becomes interesting.  The cube metaphor is of course three-dimensional, but we can also have four or five or eight-dimensional  cubes, although you don't often run into them in the wild.  So it is pretty hard.  Actually, you can find five-dimensional data cubes.  You have to do, you have to really make an effort and go to era five, whether forecast  data, these kind of things to find them.  So tables are a special case of data cubes, but they're essentially one-dimensional data  cubes.  And raster layers are also thought of, raster sort of raster data are also thought of as  data cubes because they have rows and columns and then cell values as the attribute.  Then it's the question sort of what are, what are arrays, what are the cubes and what are,  what do dimensions of data cubes?  Yeah, so I was talking about rows and columns of data, but the typical thing for spatial  temporal data cubes is that they have spatial properties and quite often that they might  be sort of aligned with some kind of spatial dimension we have, like latitude or longitude,  and have some regular, they might be ordered.  Not a property, and they might have a regular discretization, in which case we have a regular  raster.  Right?  So, so it goes on and sort of specializes and specializes.  The things I will use in these library stars that you can, the best is to install it from  get up here straight away, then you have sort of the development version.  It's also with crumb.  As I mentioned, the vector data cubes, an example is, is a panel data where we have  by state and by year we have a number of socioeconomic data.  This is a very common way of aggregated data for socioeconomic data, for health data for  instance, where you have time series over, over polygons comes off in this way.  If you would plot this, then this is said, this is the tiny layout where everything is  in columns and basically stacked and you see that, that year cycles fast and state cycles  slowest, that doesn't really matter.  If you throw that into packages like ggplot, they really like that, right?  They like the tiny layout and then they can work with that.  You can even plot a raster of this data where you have essentially here, I took the state  as the, as the y variable because then you can still read the labels.  I think if you put it along the x-axis, you cannot, you cannot read the label and put  year along the x-axis because we think often there's sort of time variability in the, in  the x direction going to the right.  If you plot this data this way using color, so we have a two-dimensional table essentially  plotted with color intensities as the, as the value for the, for the public capital variable.  And what you see here is of course, it's kind of the most boring thing.  It's kind of says, oh, there is only sort of differences between states.  Yeah.  So this is thing one and then you start thinking, yeah, well, this is public capital.  Let's look at which state is, oh, there is New York, which is very high.  Oh, there's California that is very high.  It's very unhelpful that this thing is, oh yeah, it's ordered from top to bottom to top.  But then you see it's kind of the differences between states, which is of course like, okay,  yeah, there's more people in those states.  So there is also more capital in those states.  So it's the trivial signal that you find there.  And then you're interested in basically this signal kind of overshows every other state.  First signal, yeah.  So it might be more.  What would be the temporal variability?  Are there, is there any change over time?  You see, this state goes from a little bit from dark to light, but it's very hard to  discard because most of the color range is sort of taken up by the differences in states.  So we could do that, and I'm doing that in an array approach where I convert this data  set into a stars object with the ST underscore S stars.  And this state just takes a table and then it takes it to my default.  It takes the first two columns to work out the array geometry.  So it says this is the default of this STS stars method for data frames is to date it  to first columns.  And then y decreasing is false, meaning it doesn't sort of refer to y ordering, which  would be the usual case for us today.  So I end up with all these variables that started variable three, as you see here, region,  private capital and so on.  And I move now the first two variables into a dimension table.  So this is the dimension table where we see that we have 48 states.  We have no offset in delta for this thing.  A references, I will talk about what that means.  Nothing known about point, but we have values, basically the values for this dimension are  the county, sort of the state labels.  So the state labels are preserved there.  We have another dimension, which is year, which ranges from 1 to 17, which has an offset  in delta, meaning it has a regular discretization.  So it starts at 1986, and then it has a step size of minus one.  The minus one is here, sort of it reverses the order for reasons that you would typically  do that by for raster data, which it shouldn't do here.  So something went wrong.  But anyway, it's fine.  But we now have sort of all the things laid out as a two by two table.  We could, for instance, work on that table by applying things on a certain dimension.  So what I wanted to do here is for every row that I see here is take out mean effect.  So because these states differ so much in mean, if I want to look at temporal variability,  which is my main interest, I would first take out some kind of state effect.  And I do that by dividing by the state mean.  You could also subtract the state mean, which is not our approach.  I divide it by state mean and say, OK, what I want to do is for every state, so I want  to apply for something, a function over a dimension of this array and do that over it  first.  And then I have a state and then apply functions and then the function that takes the vector  that means the time series and divides it by the mean of the time series.  And then I retain the two variable names here.  And then I get this thing out that just reversed the order of the variables, but still has  the same dimensions.  And you can see it's now entirely different values, or the values have been normalized  by their state mean value.  And if you plot that data set and also change the order of x and y and fill it with a public  capital, then you see suddenly that all of the states have the same kind of mean level,  mean color level, and that we see that certain states sort of go down and certain other states  go up.  So suddenly we see where we first saw only sort of the main effect of different states,  which is essentially the size of the state.  We now see that the temporal development by state is also varies pretty much.  So it's an entirely different mess.  So this is something you can do by saying I want to do something for all the time series,  right?  So for all the time series, I apply function.  This is functional programming and do that over a certain dimension.  And that is something you can do with apply, apply works on arrays, and sd apply works  on stars arrays.  And stars arrays are just arrays with a little bit more semantics.  For instance, for instance, it has semantics that we have a regular time series here, although  it doesn't know that these are years.  It doesn't know sort of how the time information is there.  Yes.  Yeah, always.  So in this plot, the states are different size.  And if you would like to plot that also, that you can see the proportional size of the state.  That would be possible.  Well, you see the proportion of the size of the states.  In this plot, by the variable capital.  So this is the size of the state expressed of the public capital of the state.  So of course, there are states where people are rich and people are poor.  That's not in here.  This is just the amount, the amount of something, right?  And bigger states, this is the size of the state expressed in this.  But here state is just labels, right?  So we just are interested in the differences between these.  Right.  And I was here.  So this is spatial temporal.  It's just that neither time, so we see now that we dropped the year numbers.  So neither time nor space is referenced, right?  So time was just numbers.  We need to know that these numbers refer to years here.  And space is state names, right?  So we just labels for space, right?  And ARC can deal with those kind of things very well.  But the thing is that at some stage we want to make plots and then we need to associate  these.  If we want to make plots of maps, then we need to associate these labels of the state with  the geometries, with the polygons of the states and then do that.  And then we get real spatial temporal errors, right?  So we can do that.  And we'll do that with an example of the North Carolina Sitz dataset.  So that's a very classic dataset that comes back in a lot of examples.  And had a long sort of v yet in the SPDAP package here.  And that contains data, epidemiological data.  It has several called birth number of births, then the number of certain infant deaths.  So this is the disease cases.  And then another variable, the non-white births, because there wasn't the hypothesis  that the disease might happen more or less with white or non-white persons, and was collected  for two years.  So we have essentially counties, right?  For North Carolina, we have 100 counties.  We have three variables.  We have two years.  So that is like a table with a three-way layout.  County, variable, and the time, the year.  So we can do that.  Here is the table.  If you download them and plot them as a table, I sort of used read underscore SF here so  that you get automatically get the thing printed as a table.  Some people find convenient, which I find convenient, actually.  And now I'm going to make this into a matrix.  So this is pure base R.  So I'm going to select the variables here that I'm interested in, the six variables,  because you can see that these variables are spread out over columns.  So we have the birth, certain infant deaths, and non-white births for 74, and the year  is glued to the variable name.  This is not uncommon.  We have things repeated in time, and they're essentially different columns in the table  because they were associated to the same features.  So this is, in that sense, is non-tidy because it's something that's distributed, kind of  hidden in the column name.  But that is also not very uncommon.  So we take these six variables, and then we're going to convert that into a matrix.  So that will be a matrix with six columns and with 100 counties.  But then we're going to change that in 100 by six into 100 by three by three, just by  setting the dimension of this matrix.  And now this works because the ordering of the values in the matrix doesn't have to change.  So you need to know a little bit about how this works.  And what I also do is that I set the dim names.  So that means that for each of the values of a dimension, I give labels.  So I give the county labels to the variable county, which is the first dimension.  I give this birth, certain infant death, known white birth to the second dimension.  And I give the year, as you can see here, 1974, 1997, to the third dimension.  So then I have a three-dimensional array, 100 by three by two array.  And I'm going to cover that into a star's object.  And I get this, right?  And then it sorts this out.  And it says, OK, I have a county dimension, where here the counting names.  I have this variable, which says, is this about births and infant deaths, known white  births?  You could argue this should not be a dimension, but it just threw it into a dimension because  of population numbers.  And I have two years, right?  And I have one variable called population.  So I have 100 by three by two.  So there's nothing spatial temporarily on this.  It is just a regular array with labels.  And the labels end up in the dimension values if I throw that, if I convert it into a star's  object.  So it's the way the dimension names end up is here.  But there's still nothing about references.  There's nothing about geometries.  There's nothing about knowledge that years is actually a time variable.  That will all come a little later.  What I know now is that from sort of a data cube with labels, with label dimensions, I'm  going to convert this into a spatial data cube by setting the geometry of the North  Carolina counties to the values of first dimension.  So I have this star's object here, which I had here.  And I have here the values, which are the county labels.  And I'm now replacing that with the actual geometries.  So you see now that the values are no longer county labels.  They're gone, but they are geometries.  They're abbreviated simple feature, multipolygon representations.  So the geometries are actually in there now.  So now I have an object that is spatial temporal.  What you also see is this little thing is that there's a reference system now, which  is the coordinate reference system of the coordinates.  So these coordinates only make sense if you know what they mean, right, as these degrees  latitude longitude.  If yes, to which data are they associated there, the data associated to the North American  data from 1927?  Or are they projected in what is the description of this projection?  That all goes into this reference thing, which is the first spatial data, which is the coordinate  reference system.  So it has spatial data now, and it has a spatial coordinate reference system that ends up here.  So this is something that standard arrays cannot do.  You cannot sort of go further than just dimension labels.  So textuals, it's like a little bit like what I think NumPy also have labeled dimension  or x-array and so on.  But if things are spatial, then you need to put a little bit more information in there  to make this understood and intruparable.  And these are kind of ugly things, right?  They are polygons, multipolygons, holes and everything.  So the whole terrible thing that essentially is solved by SF package is in there.  You don't have to worry about that.  So what could I do now?  Well, I would like to plot them.  What shall I plot?  I have to think.  There are 100 counties.  So I could make six maps out of that.  What I did is make three maps, but take the sum of the years.  So I say I'm going to look over the first two dimensions.  I'm going to look over the counties and over the variables and then take the sums of the  third.  Here's the sum.  Here's the function that I call as a third argument.  Apply sum to the third.  So I mean that I'm summing all the 1974 and 1979 cases and I get that up.  That means that I reuse this dimension so that I end up with a 2 by 3.  So a 2 by 3 layout.  And if I plot that, the counties are plotted for each of the three variables.  Whether this makes sense or not, but it's an exploratory plot.  And you can still see some colors, which is convenient.  You see some variability, although these disease cases are extremely low.  They're pretty close to zero in the burst.  It's kind of one thing that happens in 100,000 cases, but you still see some variability because  the default plotting scale, the color range is cut into quantiles.  So it takes essentially the whole cube, all three maps, and then it looks at whatever  it takes, like 11 different colors and takes the quantile values.  So you still see a bit of variability in the lowest map as well as in the highest map.  Although the value, you cannot read sort of what the true color breaks are here.  You would have to compute that a little bit or choose a log scale or something else.  You could do all kinds of things.  So if you want to analyze this data meaningfully, then there is again this problem that we have  states with more disease cases and then sort of counties with more disease cases and counties  with less, but then also there are counties with more people and counties with less people.  So these are just accounts, right?  So we are anyway looking, we look at the disease cases.  We are looking at a population distribution.  So you want to take that out.  So you want to sort of have the relative numbers.  So we are sort of relative to the amount of births that happen.  Where are the sort of the more certain death cases related to that?  And the usual way or one usual way to do that is to go from raw accounts to standardized  incidence rates, serve values by computing the incidence over the population and dividing  that by the mean of that ratio.  So the mean of the total population for the whole region.  So this is a job that we basically can do over, that we have to do over these three variables.  So we have to do that for every year and then basically relate the SID cases to the birth  cases.  SID is the CIs and birth is the PI's.  And that is what I did here.  So what I'm going to do now is to move, to split a dimension over an attribute.  So like a data frame, we can have, you know, we can have several variables.  In this case, we have one array, 100 by 3 by 2 array, which is population.  And it's a mix of birth counts and of disease counts.  But if we want to sort of do the ratio of the disease over the births, we can as well  split that over attributes.  And that means that we now not have one population variable, we have three variables.  These three arrays are now 100 by 2 each of them.  So 100 counties have each of them have two years.  So that's another way of looking at the data.  There's not a correct way.  There is just whatever is convenient.  And you can split and you can merge.  So you can throw sort of layout dimension over different variables as you would when  you create a data frame.  This is a list of arrays.  So it is kind of a data frame that is smart enough to understand dimensions.  And you can split them and you can merge them again.  So we split them here.  And what I did then is over-disclated data frame, I computed the sum just over this list.  So this data structure is just a list.  So I can do an s supply, which is a base r function that just looks over a list and  see what's there.  Average element of a list, I computed sum.  So I get the totals, the total birth, total sits.  This is the total birth, total sit.  So that gives me the denominator here.  The totals, sum of ci and sum of b.  So I have the sum of birth, sum of sits, total, whole data set.  And I can compute that ratio.  So that gives me the normalizing constant, so to speak.  And then I'm going to use that to compute for every county and every year.  Now I want to do this for every county and every year.  So I'm going to apply to every county and every year.  I'm going to apply, compute the ratio of ci and pi.  So this function is an anonymous function.  It just gets a length three array, right, with the birth, the sit, and the normal birth  values.  And I just want the ci and pi, so the sit value over the birth value.  So I'm going to just compute that ratio and divide it by my normalizing constant, which  was the mean ci or pi for the whole data set.  And that gives me the standardized incidence ratio for this data set.  So it's a few lines of code, and it looks very complicated, but you just do things.  You do it now for every county.  Now do things for these variables.  Now do things for every year and so on.  But keep track of that.  If you just, you know, as an exercise, do the same thing for a regular table.  You can also do that.  It's not rocket science.  It is just, in my eyes, it's much more complicated to do it the other way around.  And the advantage here is that we end up with two maps.  We just don't end up with a table that we then, again, have to link to geometries, but  we end up with two maps because the geometries are retained.  And the years are retained.  And I have CER values for the two years.  And those can then be compared.  So these two maps are then typically the input for further analysis and do questions like,  okay, is there a difference between the two years?  Is there some kind of a regional effect at the regions that are higher or regions that  are lower?  Is there some kind of clustering in this pattern?  Is there, yeah, so this is typically the input for spatial statistical and epidemiological  modeling with some kind of lattice data models.  You see here that I took a dedicated color scale.  So I took a dedicated palette, which is a diverging palette, which goes to colors from  the center.  And I took dedicated color breaks here because the value one is special, because the value  one says that I am identical to the mean.  So white cells are pretty much identical to the mean.  Red cells means that I have incidence values that are higher than the average over all  the countries in over the two years.  And blue cells means that there's cells that are lower.  So the value one is kind of identical to the global incidence ratio.  That's a normalizing constant.  So this is much more informative than focusing on maps like this.  These are the raw counts, which contain, yeah, the SIP map contains the signal of the birth  count.  So it's very, it's impossible to see whether there's higher incidence rates.  And as I said, the advantage of doing it this way is that you can work on dimensions directly  and you end up with maps that where the spatial data is still retained.  Yeah, so are there any questions related to that?  Yeah.  So, yeah.  You have the June stars.  You have the June stars.  That doesn't seem to work for you.  What do you mean?  Oh, that should actually work for this.  Oh, good question.  Yeah.  I think it should work.  Well, let me just, yeah, I think it should work because it splits the on stars sort of  looks at my roster data and my factory data and then calls the MSF in case of factory  data.  So it should, but let me sort of you, you're catching me off guard.  Obviously.  Yeah.  There's always a lot of things that don't work that might might work.  But yeah.  Right.  And then can you, you can go on and you can call a map you map you as a currency.  It's a little problem that it didn't work with Mark down very well.  So I, I, I, I'd comment, but you could look at these plots interactively on a, on a leaf  with background and see where these things are.  What the names state of the name is, if you click on a particular county.  Okay.  Let's, let's take a second dive into vector data cubes.  So, so yeah, so vector data cubes are very convenient for cases where you have not just  sort of vector data with, with attributes, but kind of structuring the attributes in  the sense that you might have a time series or you might have a set of variables that  is in some way related and you want to work with that on that as a, as a dimension.  Let's look at a more complicated or more challenging case where we have more than one space dimension  and the case there is the case that I worked out here is that of origin that destination  matrices and that case we worked out here in 432 of, of a book that is slow, very, very  slowly progressing that builds on data set that Robin Lovelace, one of the authors of  the geocomputation with our book shared.  And that is an origin destination cube from the Bristol region.  And the, the, how am I just going to, and the, the thing there is that, that origin  destination matrix basically are numbers about relationships, right?  So we have here, let's say we have 100 counties or 100 districts and we have origin destination  data.  So that means we have for each district we have the number of people moving to another  district.  Yeah, this is sort of the lifts, people live somewhere and commute to some other place  and use a particular transportation means.  We have for every transportation means we have the number from every county to every  other county.  Yeah, so we have sort of a matrix, right?  A table layout of the origin, the origin regions and the destinations and then the numbers  are numbered people.  And then we have off this table, we have several, we have the numbers where people take the  bike, the numbers for people take the train, where they take the bus or where they go  by foot.  And so that is a two, two layer, a three way array, right?  With two of the dimensions are spatial, yeah, are spatial vector.  They are the, the districts.  And here we see the district and district of one zone lighted, lighted up here with zone  33 because we will focus a little bit on that.  And if you look at the original data, the OD data, it has the origins, it's a table,  it has the origins, the destinations, it has all moving people and then sort of split  up by travel mode, right?  So we have, we have essentially the sum of these four variables here and then we have  them split up.  So we have, but they're all moving people, they're just in different modes, right?  So we have, it's, you could argue this is not tidy because it's all moving people.  We can't sort of make this dependent on mode because mode is sort of distributed over columns.  So this is one way you could organize these data.  And of course the advantage is that you only have the combinations of, you know, you don't  list probably this table doesn't list all the, all the zero entries because there might  be a lot of zeros, but lists all the positive entries.  But as you can, can guess now, we are, we're going to reorganize that.  And so we do so, we're going to do some, some, some gathering.  So we do this with tidy verse and here's the gather and spread, gather and spread are  there they'll call pivot wider and pivot longer.  What Godter does is basically that it gathers a couple of variables, sets them sort of  stacks them.  Yeah.  So it sets the foot and the car.  So we get one variable and then it repeats all the other things, the destinations and  origins.  So we get an OD with the mode, yeah, the bicycle and the number.  So I now drop the four, this, this all variable, but I'm going to stack all these four variables  into one thing.  And then I need to make a couple of manipulations to make sure that I fill a matrix or an array  in the right way, not going to go in all the details of that.  But what I end up then with a little bit of, a little bit of magic, basically, I have to  do some matching here because the, because the table is not, the table is not full.  It just might constrain, it might have dropped all the zero entries.  I'm going to make this into sort of a full array that has, that has two spatial dimensions  here origin and destination, O and D, and it has a mode.  And we see that the spatial dimensions are multi-polygons.  They're identical, of course, and they are, they're geographic coordinates related to  the WGS84 datum.  And otherwise they don't have any information.  But what you have is the values, is the mode, essentially, which is the label of the, of  the, the fourth dimension, which is the, the travel mode.  And then they can do all kinds of manipulations that are of interest here.  This is the full table.  You see that there's a lot of zeros now, but there are no NA values because of course  the NA is, will be filled with this is contingency, this frequency is the NA is will be filled  by, by zero.  And they can plot for instance here.  And what I did is do a subset and I sort of subset means the first comma means that I  am going to use all attributes.  The second means that I'm going to use all the values from the first dimension.  And the third means I'm going to use the 33rd of the, of the second dimension.  Second dimension is destination.  So I'm going now to do a slice of this three dimensional cube.  I'm going to slice through the, through the destination dimension and take destination  33.  Right.  So I get now four maps, the bicycle, foot, driver and train maps for all the people that  trend that move with this to destination 33.  And that's the highest, the destination with the highest sort of, it's a city center.  Right.  And all the people going to destination 33, where they come from, how many and what is  their mode?  What is their transportation mode?  Right.  That's in these four maps.  And you can do just do that by say, give me this plot with destination 33.  So I slice this two dimensional data, do a slice over one of the facial dimensions.  So the other one I can map, right.  You cannot sort of do the map of all the, all the combinations.  That would be, that would be a lot of maps.  Now we can also sum, for instance, this thing over the second dimensions of the destination  and then look at which is, which district is my maximum that is this, this 33.  Yeah.  That's the reason why it picks 33.  It's a city center.  And also look at, for instance, the sums of the four transportation modes, right.  Sum over the four transportation mode and get this array where the fourth dimension is  gone.  I can also do the origins, the origin totals by mode.  You see origin totals where I, where essentially summed over the first and the third dimension  and destination totals by mode.  My summing over the second and third dimension.  So you get a lot of summary statistics, essentially a very simple code by taking these origins  or destinations as a dimension on a, on a queue.  Right.  So here I took these sums and over the origins in the destination and implanted them side  by side using a log scale for the, for the variables.  So you see some kind of reason of variability expressed here.  You see the destinations.  You see that the origins are much more spread and the destinations are much less spread  and you can do some kind of for computations on this, on this kind of things.  Yeah.  So this is a, an example where you, where you have more complicated data set.  And again, you can do all these manipulations by manipulating the tables using tidy, using  the player type of things.  But it's, it's a pretty much high likelihood that it will be very clumsy and look much  more complicated to, to do it that way.  So that was a more complicated example with two, two spatial dimension, two fully spatial  dimensions, right?  So both the origin and the destinations are two dimensional spatial are, are about polypons.  So these are factor data cubes.  And before I go to roster data cubes, I want to sort of wrap up that factor data cubes basically  could be stored as data frames or tables.  The thing is that when you do with, when you work with arrays, it is immediately clear  where you have your missing values or where you have your zero frequencies to see the  account data.  So the table that, that, that Robin provided, it's not immediately clear that there are  not duplicate records in it.  And it's not immediately where the records with zero frequencies are because it's basically  it just stores the positive numbers.  Yeah.  That is that, that's the information you have.  So, so that is, that's information that's convenient to fill that basically to complete  that to fill that up.  Duplicate records is a, is a potentially a problem.  Also if you work with the same data in tabular form, then there is no index to data, right?  So the order of the records is not given.  And probably you're going to do a particular kind of analysis where the order of the records  doesn't matter.  That means there are no indexes.  If you're unfamiliar with sort of anyone who ever worked with really sort of larger data  sets where you, where you have to wait for an hour or days or something like that on  computations and it is about, it's about databases and finding data, then everyone will know  that, that having indexes in your data is, is the thing to speed computations on.  Indexes mean that you know in advance, when you pre-compute where a particular record  is.  And if your data are in an arbitrary order in a, in a very large table, then you need  to go through it and go through it and go through it, right?  This is one of the reasons why data.table, the package is so much faster because it's  pre-computed and stores indexes on variables.  Yeah.  It does basically what a database also can do and that our data friends do not do.  And if you have large, that makes a lot of sense.  Erase are just there.  The only reason they are there is that they are the indexing.  If you know that sort of, you know, if you want to know sort of the combination of the  22nd origin and then the 54th destination and then the travel mode foot, then you know  in advance exactly where to look up that value.  That's the way arrays work.  They're not just tables.  They are clever structures that where you have pre-computed where everything is.  You have that in it.  You have it.  It's a very fast way of getting to the data that you want.  And the other thing is that you can directly operate over dimensions.  You don't have to sort of go through all the records.  No, you just directly operate on the dimensions.  Say I want to sort of have the ratio of the disease cases over the total population.  And that's just two labels or something like that, right?  And also, factor data cubes arise naturally when you sample raster data cubes.  Sure, raster data cubes are up next.  Quite often you want to sort of look at what happens in this raster data at these set of  sites or what's going on for these.  How do I summarize these raster data for this set of polygons because somebody wants to  know this at a polygon level because it needs to be connected to social economic data or  something like that.  Report it.  And so factor data cubes arise naturally when you do sampling or aggregation of raster data  cubes to points or to polygons or even to lines.  Why would you not do it?  Well, tables may be more efficient.  So the data cubes that I look now at are all dense data cubes where I store all the combinations.  Right?  A clever way would be in the audience destination case that it would only store the nonzero  values, right?  And that would be a sparse array which I didn't look at.  And you can do sparse arrays, but it's just more specialized.  We often don't do that.  Things on this quite often anyway compressed.  One of the reasons not to do it is rather than be smart, use raw power.  So there is a trade off between the amount of time you're going to put into developing  a dedicated analysis that's very fast.  So the analysis would get faster and faster, but the time you spend on getting this done  gets more and more.  And so the total cost of your study is the combination of the time you need to wait for  the results and the time you had to develop the solution.  So when I talk about raw power, then you can say, well, I know how to do this with tables.  I'm going to throw it at some infrastructure that is going to do this for this big table  that's going to work.  So this is essentially the higher idea behind Google BigQuery GIS, which is a post-CIS and  steroids.  So you could throw terabytes and petabytes of data against it.  At very small compute costs and get your thing done on a table or form.  It's just a way to work.  This is also how the company Cartow essentially does everything with vector data.  They're doing everything with vector data, look at their walks and so on.  They're doing raster analysis with vector representations of their data.  They're just either the points or small polygons, whatever it doesn't care, just because it's  easier to do everything in the vector context and it doesn't matter when it is there, it's  there.  So you get things done.  So I wouldn't say that this is the future, but it is definitely something that might  be, in certain cases, the cheaper option.  Then we move on to, so any questions about vector data, cubes or my remarks about big  computations and so on?  Yeah, there is no, there is, you know, there's not a truth.  There is different paths lead to, hopefully, two identical results.  They should.  Okay, then we'll move on, yes?  Yeah, these great global grids are an interesting topic.  So I'm not really going to.  There's a person who made a package and now he's going to present the.  Richard Barnes developed a package called the DG grid, this great global grids for hexagon  grids, the covering the earth.  Together with Dewey Dunnington, I did a lot of work the last couple of months on a package  called S2, which has a different discrete global grid, the S2 grid that powers all the  Google infrastructure, which is based on a cube representation and sort of quad trees  on cube sites.  And these are essential indexes on spheres.  So when I talk about array indexes, the typical raster thing that I will now go into is like,  is it two dimensional?  Is the idea is that you have two dimensional space and have an image with rows and columns  and that that is fast.  But if you have mobile data, then the flat representation is not typically not so useful.  You could do a longitude latitude grid of a certain grid size, but then your cell sizes  along the pole are very small and they all touch each other.  So you get these problems.  So discrete global grids are efficient ways of cutting up the earth, the sphere into similarly  sized parts that can be hexagonal, H3 and so on.  And they can be nested or they could be semi-square like the S2.  Yeah, it depends on the problem.  So if you have global data, it's not a bad idea, but it really depends, I think, on how  important it is that your cells are reasonably equally sized.  If you do computations anyway and sort of, if you're able to compute the size of your  cells and take that into account, then you're also fine.  Right?  Then your cells are just un-niggly.  So it's not like forbidden to do longitude latitude grids.  It's a stupid thing, it's just very inconvenient if you're at a pole, if your region covers  the pole, it's extremely inconvenient.  So then people sort of do that and it turns the earth so that the equator goes over the  pole and that kind of strikes.  Yeah, so that's somewhat, you know, it's a very interesting discussion and irrelevant  and it's a way, there are very clever ways of indexing spheres, indexing the globe that  are good and hexagons is one of them.  You could do it this way and you just have a list of hexagons and put it in an SF structure  and it would work with stars.  It's just that stars doesn't know about the indexing very well.  It just says, oh, this is just a very large collection of hexagonal polygons.  Whatever I'm going to do with that.  But then if it would do geometrical operations on it, it uses S2, which is again, blazingly  fast, right, that has a different indexing structure.  So raster data cubes.  One more question here.  Yeah, sorry.  Can you repeat the question because of people online, people complain that they don't know  about the data.  Okay, I'll repeat the question for the back end, for the online.  Yes.  So stars and information for spatial temporal arrays.  Arrays.  Arrays.  Yes.  So, it's not necessarily spatial.  It's just in the beginning.  It's only non-spatial.  True.  Yes.  It's just we want the dimensional.  Yeah.  Yeah.  It's so, right, so repeat.  The question was stars.  What does it stand for spatial temporal arrays?  But the examples I started with are essentially non-spatial, non-temporal arrays.  They're just generic array data and that is covered as well.  It is just that it doesn't have that many benefits in that case.  Yeah, if I just have labeled arrays, I could work with raw arrays and do everything with  apply on them.  Yeah, so they then also the label, the dimension labels are being kept apply sort of propagates  the dimension labels.  So our has a rich sort of infrastructure for dealing and handling, subsetting, slicing,  analyzing, erasing, reducing arrays, permuting them, combining them.  So in the package A bind that stars also uses.  Yeah, but that and on the non-spatial, on temporal things, stars adds very little.  Yeah, there's little benefit to using it.  It just sort of passes it on to what are already does again.  So raster data.  So here is a test.drd file that I read from the raster package and I read it with the  command read underscore stars.  So read underscore stars twice to read raster data using GDAL.  And this is a vector format that actually is dedicated for our raster data sets that  GDAL can read.  And that we get here.  So we get the sort of the test.drd variable.  You see that I get two dimensions here, x and y.  So I now have two spatial dimensions, yeah, two dimensions, x and y that kind of together  sort of describe space.  So this is a very extremely different situation.  I have x and y and I need both of them to sort of find a certain location.  So in vector data tubes, I had just a sequence of either points or polygons or lines, whatever,  that each of them described two or three dimensional space.  Here I need two dimensions to essentially to get that.  I could select a value sort of for x and take x with value three and then I would get a  constant x meaning I would get a column of my raster out of this.  I could select a column or I could get a column subsection of that raster by selecting x values.  So we see here that several things.  We see that x ranges from 1 to 80.  So it's an 80 by 150 array, but that we also have offset and delta values.  And offset means that is where the raster starts on the x side.  So the x edge of the raster is the offset and the y edge is the y offset.  And since y values, y delta, so delta specifies the step size, so the grid cell size in x  and the y dimension, you see negative y values, negative cell sizes in the y dimension for  the reason that the y coordinate increases when the row index goes down.  So they have opposite directions.  So that means a negative cell size.  Why goes up going north?  The typical images have cells that the image, when you plot the image, it goes from top  to bottom.  This is the way image data is almost always stored.  You don't have to do that.  It can also be the other way around, but almost everything that is stored in file is stored  that way because PNG, JPEG, everything does that.  Images start bottom to top.  So that is the delta.  And we have a reference system, which you can see here is this 4 string, which are now  deprecated.  And the other things are missing.  We have a special label that says, is this one of the x and y?  Yes, this is the x and this is the y.  The reason is that not always x and y are called x and y.  They might be called latitude and longitude or lat and long and so on.  And somehow they meet, they're mapped to what is going to be plotted as the x and what's  going to be plotted as the y.  This sounds trivial and it's not always trivial.  There are a lot of cases where if you have a polar projection, it's not clear what is  worth it.  And then if you plot it, we use this theory, this scale, and we use equal breaks and we  get the usual plot of this new data set.  So we see here a couple of things.  Raster files are data cubes.  The attribute has the values for all combinations.  So also the white cells here are filled with NA values.  You see here that I have 6000 NA values.  So this is all sort of the whole roster is filled and I just plot, I just show the colors  for the non-NA values.  Their vector data cubes have a single spatial dimensions.  Raster data cubes have two spatial dimensions called with x and y labels.  Or they might have three, if it's a three-dimensional data cube with foxholes or some kind of meteorological  data and so on.  This case is an example of a regular grid.  That means that offset and delta are filled and that indicates that we have constant cell  sizes in each of the directions.  Here there are also square cells, but as you see they don't have to be square cells.  This one might be 40, this one might be minus 50.  It's fine, but they're constant over the space.  That means that I can compute my coordinates from the index by this equation.  If I have offset and delta and I want to know the coordinate of index i, then I can compute  it this way when my index is one-based and then for the value one, I get the cell edge  and for the value 1.5, I get the cell center.  We have the cell space, which is continuous, and the index, which is also then thought  of as continuous.  This took a while before I understood.  Then why is the delta for y is negative?  That means that the coordinate values decrease when my index increases, when I'm going down.  Also, as I mentioned, you have a reference system.  That means that I can link this data to other spatial data and put it in some place in the  earth.  What is plotted are the non-missing value cells.  What it can do with this data is convert them into an SF object.  That means make them vector data.  What then happens is that, and I say, don't remove missing value, which would by default  happen just to show you.  I get polygons.  I get here then 9200, which will be the product of 80 and 115.  Polygons, and if I plot these, then I get a plot like this.  You see the little polygons plotted, where I also plotted the polygon borders in gray.  You see all the grid cells plotted as small polygons with a gray border, and the values  filled whenever there are non-many values.  You get the same map in a vector representation.  You can equally do that.  This is what I meant by what Carto does and what BigQuery GIS would do.  Look at the whole thing as a bunch of polygon data.  It works also.  Of course, you wouldn't do this with massive grids.  This works nicely because it's relatively small.  If you're going to do this with very large grids, then you're going to set your computer  on fire.  Not literally, but it will require so much memory that your process is breakdown.  These things actually become really interesting when I use multiple layers, raster data cube,  so-called raster stacks.  Yes, your question?  Just a brief question to something.  You said x i, where you indicated x i with the offset grid.  So that's the net edge.  There is a value offset and delta that you find in this table, offset delta.  You filled it in in this equation, and you give it an index.  Sure.  I know about that.  That's the net edge of a box, say, area.  And why?  Is it the top reference of the box?  The question, repeat the question.  For x, I can see how this works.  It's the left side.  For y, is this the top or the bottom?  And here is the magic.  It depends on delta.  Delta is either positive.  If delta is positive, then an increase in index means you get an increase in y.  You go up.  And the offset will be the bottom edge.  If the delta is negative, taking an index of one, the first thing, cancels this element,  means that the offset is the top edge.  So the meaning of offset depends on the sign of delta.  That is the trick, how to think of offset and delta.  Thank you.  Good question.  Right, multiple layers.  So here's an example from the RGDAL package, which has a familiar JPEG image.  You can also read JPEG images.  Reading JPEG images is not a problem.  You get three bands, RGB bands, because it's a JPEG.  And you see that here are dimensions.  And you see that these are regular discretized.  I again get an offset, which is 175, and there's a delta of minus one.  So again, it's going down, indexing.  But the reference system is missing.  So it's regular, but this is not Earth data.  This is just an image.  And that is also why the cell size is one.  But JPEG can read such things, nevertheless.  But it's just unruns you in spatially referenced.  You get the intensities in values from 0 to 255, which is typical for JPEGs for the three  bands, R, G, and B. And then you can use the function SD and RGB that computes RGB values.  So R has values for colors.  You can use the RGB function to compute a, ah, it needs to be maxColor is 255.  It gives you labels.  This is text data.  This is character data.  It gives a hex representation of a color that can also include an alpha value.  And this basically throws these three bands together, or possibly an alpha band as well,  and converts that into color values.  So that reduces the third dimension of the band.  And I get these color values out.  And when I plot that, you see the regular RGB thing.  I don't know why it's not plotted.  It's the thing you've seen before.  Let me try.  Here it is.  So I need to.  I'm not so good with.  I'm not so good with our studio, to be honest.  This is very embarrassing.  Oh, yeah, this is it.  Yeah, thank you.  Yeah.  Thanks.  And then you get the unsurprising, you get the color map of the R logo.  And that puts RGB on top on computer colors.  Right.  So this is something you can do to compute colors from values.  Right.  Multi-band data is more common with your handle satellite imagery.  Here's an example of an excerpt of a small Lanzon 7 image used to cover part of Brazil.  You see it has six bands.  It has all the 30 meter bands.  And they're there in UTM zone 25.  Yeah, so Lanzon data is in UTM.  You also see here that the point label is set to false and no longer to N a.  That means that the data set that we have been reading here that this stiff file has metadata.  That tells you that pixels are actually the pixel values are associated with special areas and you're not associated to point values.  Yeah, there's a difference between one and the other.  You could have other data sets like, like the M data digital elevation model data where your pixel values refer to an elevation at the cell center.  That is a different that has a different semantics. If you then would sort of zoom in, you would like to interpolate values rather than draw them as square boxes.  Yes, here you hear interpolation is sort of harder to imagine.  Of course, this doesn't mean that the cell values are a constant value throughout the cell.  They're just associated to the square and they're of course they're confolutions because this isn't this is an image sensor that that sort of measures averages over small areas that are not even square.  Right. And if I plot that, then you get you get a plot like this, and it essentially histogram stretches over these seven bands, right.  So it takes all the values of these of these array. So this is a 350 by 350 by six array takes all the values, and then takes the quantiles and sort of assigns course that means that for these six maps, the sort of the 10 different grade  scales are kind of equally distributed. Yes, so it's similar to what you do with high dynamic range photography that you kind of distribute colors equally not get very large, light or dark areas that you see something if you would do it in in standard constant or equal  interval range the way that roster does it and most of the times all all the maps will have the same value will be all dark or light you don't see anything. Yeah, so you can override that here I did the same thing but then I didn't do this with a joint Z limit so  here all the maps sort of refer to a single color scale but I can also do the histogram stretching for each of the bands. Yeah, so I do that here so here you see more color stretching and of course no key anymore because that each subplot now has a has its own key  and has its own color breaks. So you can do this to sort of blow up variability. And it also maximizes spacing of bottom region so you did you see something right instead of putting legends everywhere.  Right, so you can then also do RGB composites here this is the true color composite this is some kind of false color composite by by computing by computing RGB this is by plot and using argument RGB I could also do as the RGB with the band selection, and then  compute the RGB was different ways of doing this. Transforming and warping rusters. There is the thing that if I want to sort of combine Sentinel to and then some eight and it might be the case that these two don't have the same code and reference system is definitely the case if  you put something on top of modus data which is this sign is vital grids, then you need to work with that. You need to do something and this to illustrate that I created a little data set here.  This is this bounding box. Yeah, this is pretty much sort of a large section of Europe. And then I said, converted into a stars object, which has essentially an empty at has zero values. And then is a 30 by 20 grid with with one degree grid cells associated to WTS for covering  that part of Europe. And then I did this not going to talk about what that means. Maybe in break and loaded some our natural Earth data, which is which is global Earth data with some, some, some, some parameters, some parameters and downloaded the Europe, the Europe  section, which includes Russia, and then computed population densities by basically taking the population estimate and dividing that over the area of each of the of each of the polygons this is polygon data.  And then I convert the units to one over squared kilometers yet to not get very crazy values, very crazy low values because the this is one per square meter.  So I put all these data, and I'm next going to rasterize that that means that I'm going to take the population density polygon data and rasterize this over this over this object, it's that they just made which is the one by one degree grid. And it just sort of takes the does the  population density values for these grid cells, right. And I do that. So it's very course map of a long acute latitude. You see that, although the grid cells are square in terms of degrees.  Yeah, you think square grid cells one degree north and one degree east. Of course, one degree north and one degree east doesn't everywhere on the world mean the same distance. Yeah, so you see now that the plot methods, given geographical coordinates taken aspect ratio so that in the  map roughly one meter north equals one meter east, which is not the case with degrees. Right so it doesn't aspect ratio correction on as what is essentially an equidistant equidistant rectilinear projection.  So it's not degrees let it as if it was x and y, but an aspect ratio correction to make sure that countries have the kind of the shape that you're familiar to and they're not like like Germany is this compressed term this kind of round Germany which you see when people ignore  it. So, so this is something you do and if you do that then you also see that the roster is no longer has no longer square grid cells. Yeah, but has grid cells that in kilometers in the center of the map, roughly, roughly has one degree in one  degree and see that one degree latitude is like point eight degrees long.  That should be the other way around. Okay, so what I can do what I wanted to make the point now is if I if I'm going to transform this to another coordinate references for instance the ETS 89 Lumber equal area projection which is the  recommended European sort of plotting which is at least equal area plots of Europe, then I can do this for this grid for these great values and I get something that is that now looks the same is the same data and is still a 30 by 20 grid, but has a different reference  system and has a big bunch of coordinate values and it's now it says I'm a curve linear grid. Right so what I essentially did here was nothing but just change the coordinates of every grid cell.  Now here I knew the coordinates because I could compute them from my regular discretization. Now I'm moving them to a new system which is either where I have to compute for every grid cell, I have to compute the coordinates in the new coordinate system, and they look like this.  So it's still a regular 20 by 30 grid, but every grid cell now has its coordinates in this in this coordinate reference system which is EPSG 3035 which is the ETS 89 LAAEAEAEAEA projection that actually every European citizen should be aware of that this is the way to plot your European to use for your European maps.  So it's just a little bit of a digital familiar. All European publications will be will use this. I think.  Right so it's still a grid. It's just not. It's no longer regular. It's no longer regular. It's so called curvilinear grid. And that is that can be convenient or inconvenient convenient thing is that all the grid value are still there.  I just changed the geometry. Yeah, I just the geometrical transformation. And of course, at some stage you want to work then with our systems that need regular grids and so on.  Otherwise you can then do is warp this grid warping means that I have some kind of up shaped thing and I'm going to burn it into a new regular grid. Yeah, or you can work from anything to anything. But this is what's happening here. I take this curvilinear grid. I'm going to make a target grid.  This just the sort of the default grid sizes that takes the bounding box of my population data on this reference system. And then happens is that it's going to redo the values. Yeah, so here you see then the plot of that new grid that you see here the sites, you see the ruster cells and much much much smaller.  And of course in the values you still see the lucky thing of of of what I had here. It's just rusterized to a new roster which is a regular roster in ETV RS 89 L a EA. Yes.  I don't know how to open the chat.  Yeah, then I get back to my thing.  Yes, you can load the news data set. That's a good question. So, if you don't have this data set you may have to install the roster package. Yeah, because that is where it's going to read this from this data from if you install the roster package you will not have no problems to do this.  Okay, good question.  So this is what warping does. Yeah, so it works one sort of roster in one corner reference system to our corner preferences.  And it means it's going to recompute values. Yeah, so the totals of your roster will no longer be the totals. So for density, they still makes some sense for population density. If I had done this for population numbers, it would have been entirely senseless.  So it should be a population here over roster cells, and then re rasterizing that will sort of destroy the whole idea, because they are kind of extensive variables. So, in that case, you would use as the SFST inter interpolate area way that use area weighted interpolation to go one from one,  from one rasterization to another one and do that on a polygon basis. Yeah, the way Carter would do it.  So now it, yeah, sorry, any questions to related to the multi layer rosters, curvilinear grids grid warping RGB.  Then we will go to the more exciting things, which is roster data cubes, time and multiple attributes.  So we haven't seen any time variability here. Now we're going to load a data set that is small data set it's available in the package stars. That's a net CDF data set right net CDF is a format that has been around.  It was, I think it was one of the first sort of libraries that I discovered online. This is in the, in the, you know, decades before HTTP was invented is in the time that we used ftp to, to sort of discover data sets right you would ftp to some first of all you would go to,  to groups to use net groups and see what's there right to read emails and sort of mailing lists distribution system and see where you could find things and there would be no outsides or something like that.  And they would have net CDF files and a library called net CDF that way you could read and write these things right this thing still exists and has sort of developed over the years and is still dominantly used by modeling communities right so it is.  They use it for spatial temple they use the grid of data.  And, and for data where you have multi dimensional grid typically time varying things like weather data, like Earth observation data a lot of things are nets.  So here, this thing is being read and you can see that it reads to sub data sets or to sort of attributes in this data set, and it is a three dimensional grid layout, where we have an X and Y.  And you see that the deltas are one eighth of a of a of a degree.  But you don't see a reference system. Yeah, so that's unclear this this is degrees, let it be long queue but it doesn't say so.  And then we see that there are two value precipitation and something with temperature.  And you see also that that it reads units. Yeah, so these data sets come with actual measurement units, which is which is very cool.  And they are being read and they're being propagated using the units package, which is another story, which, you know, which is extremely, which is an extremely good idea.  But occasionally causes some headaches for people because you sometimes you want to, you need to get rid of the units before you can do a next step. Yeah, if you want to put things with Gigi plot that not always works directly.  Here we plot this data, and you see what it plots is a sequence of maps that loops over the last dimension. Yeah, so these are 81 by 33 roster data sets, and they have 12 of them and you'll apply them in the labels are the dates that are associated with each of the  sets. And these are apparently monthly something monthly averages or monthly totals of precipitation, probably total precipitation is an average temperature.  Otherwise, it would make sense. But we would have to look in the metadata of this data set to get more information about who created it for what purpose what it is exactly I think this is climate down scale climate data.  Another example here is where we have data that are distributed over that over several net CDF files. And this is typically by monitoring programs where a satellite flies over and measures something and somebody  realizes that data and sort of sort of converts it into some kind of product like sea surface temperature and uploads a new file every day. Right. Instead of sort of rewriting and very large pilots and larger and larger, the file is added every day and you see  here that the that the date is part of the file name so that makes it easy to to find these things. Yes.  Oh, thank you. Oh, that's good. Yeah.  Thanks a lot. It should work.  Yeah, it's half an hour to go.  Other questions or warnings.  I have a set of net CDF files here and it's going to read them in a single bunch and it's it the read stars here is going to sort out that all these are also data sets that seem to be arranged along a time axis it's going to go into each of the file read the time stamp for it and then sort of collate  into a regular time dimension where you see that I have nine time slices and that this is the first value so this is an offset this is a post CT structure so you see the post CT value of the of the first day which is two hours Central European standard time so it's midnight  UTC just my system wants to plot it in in my local time zone, which you can change and has a delta of one days. Yeah, so it sorts out that this is regular time this is a regular time axis.  And it doesn't know whether this is associated to time point values or time intervals that is not in the metadata. Although net CDF actually can have that as metadata it can if you do temporal aggregation, you can say my temporal there does time bounce and sort of it's really an average or a maximum over a certain time period  and then it is weaker on the on the spatial thing because it doesn't have a current reference system for whatever reason that is probably encoded in a way that is not interoperable or that she doesn't doesn't recognize.  That doesn't, she doesn't recognize this as degrees in some kind of associated to some kind of date. What you also see is there's a fourth dimension which is called Z level. Well, this is C surface temperature.  This is a surface right so it, they're not much interesting going on the elevation zero, it is C surface. So this is zero meter. So it has essentially a unit for the elevation, but it's a single layer elevation.  So we know it has metadata. This is on a certain elevation that I know zero meters with respect to C level, but it doesn't have any variability here because it is surfaced.  It causes a bit of trouble sometimes you want to plot it. It's going to loop over the third dimension and so show, see, show the sea level which is not terribly interesting.  More interesting is the variability of time. So I can drop this. I can use a drop to drop singular dimensions dimensions that have only one value.  It's a, it's a 1440 by 720 by one by nine array. But if I dropped it one dimension that has only one value, it remains the same information.  I just lose the information that this is at zero level. So I have to, so it's there it's in the data. It's a four dimensional array, but with so to speak and obsolete third dimension.  I drop drops that one and retains the other three that have variability. And then I could plot it. And here is an example where I used Gigi plot and GM stars with a couple of other things to to plot these data and the various thing and so on.  So basically you, you give it, you give it Z and then see one means the first attribute. If I don't do that, it takes the first attribute anyway I have a certain alpha value, and I'm going to down sample this yeah go to down some of that means that, although I have quarter  degree, quarter degree grid, I'm going to down sample this to a two and a half degree grid, because my plot is anyway so small that every single pixel is probably already five degrees or something like that.  So, I don't have to do that. But I'm doing that because otherwise it takes the eternity before this plot is finished. Otherwise it's going to do all the pixels.  Yeah, and then it takes like five minutes and every time I compile my slides to wait to five minutes. I don't know if it's five minutes but it's it's a factor hundreds that I'm reducing here in the number of pixel data show.  And in the end, the thing looks the same. Pretty much the down sampling is something that we haven't seen before, but it also happens here it happens in every plot.  So, I have a lot of different commands, but these are base plot commands and there I can, you know, when I do a base plot can check what it was my device how many pixels do I have how many pixels is area have here I mean ggplot and you can no longer do this because you are sort of lost in the ggplot abstractions of everything.  And at the moment that ggplot, you know, ggplot distributes tasks at the moment that you're in a geomunders core something you cannot sort of ask ggplot how many pixels does my device actually have.  You cannot do that. Yeah, interestingly, I had, I had discussions with some of the real sort of ggplot gurus.  Thomas and Claus Wilka Thomas better than so on and you can't you can't do this.  And it's very, which is very regrettable that you can't do this. Yeah, so sometimes ggplot will take eternity to finish. Yes, and then remember that you can do this down sampling by hand, regrettably, which you could do automatically if you knew if if if you know that.  So this is developer talk, but also warnings that some things, you know, are difficult to do in our other operations on data cubes are subsetting we already saw if I take the first sort of a single index.  Then what happens is I just do a listing. Starsh objects are lists with a class in a dimension table in an attribute. There are lists of arrays. So if I do a selection of the first element, I just get the first list of the second list element which is here, the anomaly.  I just drop all the other arrays and keeps the dimension table because that doesn't change. Other things is that I can select other elements by saying comma comma means that I select all the first attributes.  I don't specify anything I select everything. So it selects attributes is the first square bracket selector, then it selects first dimension second and so on. So first and second dimension x and y. So now I selected a subset that has the first 10 rows in the first 12 columns.  And then it's reversed the first 10 values of the x variable meaning the first 10 columns, and then the first 12 of the second dimension was the first 12 rows columns rows. This is how the thing is organized.  So there is this, there is this fortune about illogical, which is, which is one from Brian Ripley in 2004.  The medical conventions are just that conventions they differ by field of mathematics. Don't ask us why matrix rows are numbered down but graphs are numbered up the y axis, nor why x comes before why but row before column. Right. This is what I just, what I just mixed up.  Yeah, and you always a row column, and you say XY, but it's calling wise row. Right. So if you talk to me this, this is just conventions.  So I see level is still in there and retainer. So, and you can also select and see level of course they only have one there, and I can also select time.  Yeah, that means that I need four commas right for the attributes x y and z level, and then I get a smaller time slice right I select from time time item five three from three to five.  And now you see why the from is there, because I selected from three to five, but I still have the same offset.  Yes, and if I would say okay, those are three time instances, I do from one to three, then I need to re compute the offset. Yeah, because the offset shifts.  And that's why this is why this from field and this is sometimes confusing there isn't st underscore normalized that that does that. Yeah, that will sort of change this into from one to three, and then modify the offset value with the appropriate  number of deltas.  So, so it can also be used to shortcut intersections as like here, if I do a global grid, and then I plot sort of this, the section that intersects with that intersects with the Europe, European polygons.  Let's see if we can do that.  I think I didn't do that because it looks very ugly. And it looks very ugly, because, because any is is European data set that includes entire Russia and that crosses the deadline. So you end up with a grid anyway that is global because the deadline is involved and all the kind of  I didn't actually, I didn't actually do this.  But the idea is there right so you have you do an intersection you select by, by sort of the area you basically crop the area by a certain geographic region that can be a simple feature geometry.  So you do is cropping or making sort of a sub sub roster with smaller area, or you crop by a polygon set, and you then basically and burning and a values in all the cells that are outside this area.  So you have the verbs that tiny study verse type verbs that work with stars object. There is a slice that selects dimension values index. There is filter that selects dimension value using the dimension value so you filter you use coordinates for slice you use index you use  there was pools that pulls out an attribute select that selects attributes and mutated computes new attributes based on existing ones. Yeah, so you could mutate and then say compute sst times five or something like that or some conversion into Calvin or some conversion  into far and height or something like that do that on an area based basis right for all of the cells. Then there are also kind of these kind of map algebra type things we can  load in the  set from from the from the launch out imagery, and again, for instance, multiply by two, and it gets something that go to 255 because to 510, and it can also do things like mask things out so so for every for every cell that is smaller than 125 give  you a value in a or something like that, and then do is basically do these kind of array array logic right so you basically retain all the cells that have values large under 25. So a number of sort of things that that that you could do that you think  you could do with these types of things like algebra like, like expressions you could actually do with these stars.  So I just saw the sdssf, which basically takes rasters and then draws either points on all the roster cell centers, or draws little square polygons around all the rosters, depending what you want. That is one way.  There is the raster eyes that we saw that goes from vector to roster and raster eyes is then that essentially calls the GDAL utility rust. It's not. It's part of the GDAL library, the raster eyes routines there.  And then there's the issue of proxy objects. Yeah, so and so there's a couple of ways of going from raster to vector and vice versa.  There's also aggregate cells over polygons and computes whatever function to some or the mean over that polygon of all the values.  There's even a very fancy one that does that using the package exact extractor which is a terribly fast and good package for for computing aggregates over polygons that basically looks that cuts all the little square raster cells into parts and sort of use the  part of the raster cell that falls into the polygon part of false hours. So you get this really, really area weighted aggregations for some some means. And that does is an incredibly fast way with a very clever algorithm.  Exact extractor.  Exact extractor by Daniel best one of the SF contributors and Geos developers.  Then there's the thing how do I how do I work with large raster data. Yeah, so I have my hard drive is full with Sentinel files and so I'm going to go through that package raster does that and in in a particular way that is not not too transparent.  We do that in stars with do that in somewhat different way in sense that we have this idea of proxy objects. And I'm doing that here with with Sentinel to  say that you can use a granular or tile or seen or whatever you want to call it that is available in the star state of package. So I'm going to create here this file name which is a very long file name which has a prefix, essentially to tell GDAL this is a Sentinel to L one C data set and tell  you where it where it is and then give sort of the long thing of it need the thing that it needs to read inside that zip file to get a metadata, and then it reads it. Yeah, so this is all, you know, this is all something you need to do once and then for every  time it works the same thing. So it's going to read it here. And what you see is that when you get back is the stars proxy object with one attribute in this file. And this has the four 10 meter bands. You see here that bands are named before the five p six and the eight I think  is the one that I have 11,000 by 11,000 pixels, each of 10 meter size, and they're in a particular UTM zone.  So this is the four 10 meter bands in a Sentinel to tile. And it's not reading anything. So this takes like less than a second because it doesn't read anything. It just says the data are here. I know where the data are. And this is the metadata dimensions. Yeah, so it tells me sort of the dimension of this area without reading it  because reading it would would sort of take a lot of take most of your memory. We can control whether this is being read as a star's proxy by setting proxy is true but at current state of affairs. It's being read as a proxy object by default because it's so large.  So it's first going to read the dimension of the thing. And then besides, okay, there are more than 10 to the power eight cells. I'm not going to read this into memory because most likely you don't want that.  I'm going to take the have this as a proxy object. And then it tries to do lazy evaluation that means that it doesn't do anything up to the moment that you wanted to do something. Right. So at some stage you want to do so you so you can compute for instance here.  And the VIs by applying a function called any vi that called that computes your any vi values over the band dimension to the accent to look that over the X and Y dimensions right so apply it to the bands.  This is this is the expression so band four is near and then one is read.  So if you do this, it's again it's done within, you know, within a microsecond because still it doesn't do anything because you didn't ask for pixels. It just says okay I still have this property object.  It just now has a call list I need to do something. So if you want to have something and going to do something. Now if you plot this thing. It also does that in like something like seconds or so.  Because because what it does is basically it computes as many pixels as you have computers many cells as if you have pixels on your screen or on your device.  Right so you can do this on a 10 meter resolution computer and the vi's and it would take like 10 minutes to read the whole thing and then compute it and probably write it back because it's too large to hold the memory and so on.  So it would take 10 minutes but if you only want to plot this then you basically look at the plot of what will it be 600 by 600 pixels or something like that.  And I'm now going to loop over the 600 by 600 pixels read them compute any vi's and plot them.  Right so this is a bit of a trick so so this plot that any vi actually looks at this object and say okay there isn't this object but it still has homework to do, but I need to plot.  So let's look at what the homework does. Oh it's an sd apply or it doesn't do anything special it only works on the bands so I can do it for just as many pixels as I need is no filtering going on or if it was filtering I of course couldn't do this trick.  Yeah so it looks and says okay I have to do something and do this at the screen resolution and then plot it and that speeds up things from, as I said from something like 10 minutes when you would do it for the entire day to something that this within seconds done.  And the trick it does of course is that it uses it in which permits is usually over views this is rush to tile comes with lower resolution overviews of all of the bands. And if you asking sort of a down sampling of whatever effect of 600 and it can use your work on one of these  overviews and only read a very small amount of the actual data set. This is a pattern that everyone who works Google Earth engine is familiar with. Basically, you put it, but that does is that it just computes the screen pixels.  And then you zoom in, it's going to read compute the screen pixels. So although you work on massive imagery, you just do computation that are sort of proportional to to the plot size you want.  And finally, a more complex, a more real example that that I sort of cooked up in the last few days because I discovered that Mario's apples material for tomorrow was was online.  So I know I'm always curious like what's going on there. So my years will do a sort of a serious sort of quasi serious scaled data cubing.  And it's not going to work with petabytes, because you don't have petabytes on your on your laptop, but it's going to work with a sets the size that you can still download it and work locally to get the idea.  And Mario's compiled a number of lots of eight scenes over for a particular area over a particular time range and creates cubes from that.  And I'm going to give a very poor man's poor man's version of that. Yeah, so what I do, for instance, is I look at only at a single scene so I can I don't have to, I don't have to worry about all different places and times and so on that typically satellite  scenes they're not they're not in the cube right so it's an orbiting satellite so it's a trajectory of images you can. The stars hasn't as the mosaic so you can mosaic things together but you then need to first sort out which things to mosaic.  And that is exactly what G dot cubes does what I sort of did not do here. But then I downloaded Mario's data set, which is, which is here you can you saw all these all these tiles here which are in directories and  in in one of these directories you see sort of you see a number of tiff files. I'm going to. And we're going to load these seven bands of these of these duties.  And do this in this particular scene. So, I'm going to match all the file names that that match to this patterns or band one to seven from essentially from this, the directory that starts with this expression.  And then you read all this so that gives me like 120 files. I'm going to read all of them in memory so not as proxies and memory. And then I'm going to warp them I had to do that because there was tiny differences in there couldn't sort of put them as a pile on top of each other.  And then I'm going to do a slight differences in in the number of rows and columns so I kind of want them anyway to to the geometry of the first one right hopefully not too much lost.  And then I'm going to glue them together which can be done by the sea basically I do a sea with with the 120 objects, combine them and and then I get essentially that that set.  And then I got the time series actually from the scenes somehow right because the time. Yes, you can see here, the time is encoded in the scene name.  So this is going to happen here I'm extracting this, this little section of the file name, and I'm going to read that with this format string, and I get my times.  Sort of obscure way of doing that. I wouldn't recommend it to anyone, but it's possible that then I'm going to take the band name from the, from the following that takes this little string.  Right that tells me the band name. So I get band 127 and so on. And I'm going to merge this whole thing. And then I see that I have a seven 60 by 70 70 roster with 182 layers.  And then I have to assign band and time because this is these are the two layers are seven bands for like something like like 40 different time instance.  So I need to then read and mention this thing. And what I sort of do and now comes the ugly part. So, so this is why, you know why you need proper packages that do this in a sort of in a readable way.  I'm going to sort of just work on this raw dimension table, remove this attribute and then create two new dimensions.  Actually, one with the unique dates and one with the unique bands, and I'm going to read and mention this thing. And then you see the dimension table looks like this.  So I have now 26 dates. Yeah, which has all date which apparently have no offset and delta so they're not regular. They're not found to be regularly distributed time.  And you have seven different bands. And the question is whether everything is correct now and I'm not entirely sure whether everything is correct. I just did this.  No, everything is correct.  This might look okay so here I do a slice over date so I take one particular date. So the first date which will be the 6th of August, and then you get a plot of the seven different bands for that date so this time slice.  Now this was the essential him thanks by that they did all values more than zero were missing. So I had to assign them an a this this you can do this for a four dimensional area like this.  And here I did a plot of a time series of the first band. So this is band one and a time series and I'm entirely unsure whether this makes sense.  But anyway, you can you can get an idea time slices, then slices, and look at these data as a data cube and sort of get through them.  And you can even compute NDV eyes by doing the same trick with them with then five and four which apparently correspond to and near infrared and red and apply this to the first three dimensions x, y and time.  So we compute any vi over the we use any vi function to reuse the band dimension and we get back three dimensional data cube.  It is not just time series where where band was reduced to end the vi. And we get these in the vi time series, not sure whether they really make sense.  But this is what what came out. So it's still kind of work in progress. But just to show how you can manipulate these sort of.  This is all in memory. But in principle, you could do all these kinds of things out of memory, but and sort of still have plots exploratory plots and do simple functions on dimensions.  In a relatively responsive fashion. So without doing the full computation on the on the entire.  I hide dimension data cube, right, you see very small surplus. So the number of pixels that are being computed is here is also very limited to just kind of the same number of pixels that you see here maybe 100 by 100, rather than the rather than sort of the 700 by  700. So that is a kind of a factor 50 speed up that you get in the plotting, plotting action.  So these were some kind of, you know, some demonstrations of how you can juggle around with with factor data tubes with raster data cubes, where it really becomes exciting if you have more than four dimensions, more than three dimensions.  So more than three dimensions is also the, the use case where, at least as I understand it, the roster package and its successor.  And kind of where you get stuck, because it's, it has a roster stack, as it's, as it's multi dimensional thing and if you have a stack of bands and times, and you need to do the whole administration like where are my times where are my bands, and you can you can do a function through that.  It's not two dimensions, right. And the data structure you can handle is really three dimensional doesn't sort of work at more dimensions than that.  So that was one of the motivations to do this work.  So are there any questions from your side, which I can take and then repeat so nobody in the field will hear it.  Yeah, for example, very efficient overlay for the hard to do it.  Well, as I said, if you do, if you have a complex data cube object, like, like here, like in this last example.  So, when I do this, when I do this plot, so if I, if I use stars proxy object so here.  At this moment everything will be computed because the thing is in memory, right the whole thing fits in memory it is in memory.  But if it starts proxy object, then it has this lazy affiliation so the object is on this might be massive, but it looks at the order of operations and says okay you want to plot plots how large my plotting device how many sub maps and I'm going to plot how many pixels  do I need for a single subplot. Oh, that's 100 by 100.  So I'm going to sub sample and just read at a lower resolution. Yeah, and that goes to GDAL and GDAL is very fast and doing lower resolution operations and you don't have to tell it.  It does this because it knows it's a plot I only need so many pixels.  And does the star the same thing.  I also don't know there's very little information about what they're actually doesn't really doesn't do.  I don't know what overlay means.  Extract.  Well, you know, there's a lot of things here are also done in GDAL which is C++ and that's what they're also does does calls these C++ calls C to retains.  I hope it does, because we understand what they do.  But it also does a lot of things in ways that are sort of hard to track down the exact details.  Well, this is not I mean it's just unfriendly to say I'm very sort of, I'm a big admirer of Robert Highlands works. We looked at RGB composites of plots, and I looked at what the map view code did.  To do this, and I looked at some other package that also did that, and they all would sort of look at RGB values and then take this, this two and the 98 quantile and throw everything outside of that away, basically set it to 289 quantile and then use that range to compute RGB values to do a stretching.  And that happens to come out of raster code, which is not documented way you can't control the quantiles, but it does that.  And it doesn't show up in the help page of the plot RGB commands. So plot RGB does that and if you want to reproduce with plot RGB does you can look in the code. You see that it does that.  It doesn't explain that it does that and doesn't explain why it does that. So these are the kind of the little things that that kind of start to live their own life and then people copy that because they think what RGB is the thing that is true, right, because that was their first and everyone uses that.  But I don't know why it picks the second and 98 quantile and not the five and 95 or the one in 99 or the whole range. Yeah, there are always like funny values at the tails in a lot of satellite imagery, funny high value so if you take the data range as a default which might also make sense then you get completely black images.  And so there is there is stuff that you ran into which is kind of came up for that probably Robert got his inspiration from somebody else, but we don't know and it's not documented.  At least they couldn't find it.  So that's the essence.  I don't, I don't know what all the functionality of the roster packages.  That is those are two common use cases that is what it does. Yes.  So if you look at sort of the, if you look at the amount of functions that the roster package provides it's, it's massive.  A lot of multi core things. I think the, I think that as the apply has a couple of options, because of lies, always like in independent job because you work on dimensions.  The chunks are always independent. So, as the apply has had options where you can pass a cluster and tell whether to use future or to use water parallel apply or something like that.  So, I guess parallelization options.  So, you know, you can use the roster object calling vocal on the roster object and converting back so that was like a three liner or so that you could throw in sd apply because the focal is essentially over x and y dimension, you leave everything else intact.  So, for each of the time for each of the bands you do the focal on on the XY slice. Yeah. So if you put, I think I'm, I think I'm offline, because I did some VPN stuff that they shouldn't have done.  So if I look at stars issues and then search for focal, focal filled with stars.  Here's an example.  So you do an sd apply with with the focal function. Yeah, that doesn't.  So, the focal and roster. So does the roster conversion the fly and cost.  Yeah, we could port this code is this is also kind of kind of the same thing the question is whether how this will work with with massive with massive grades.  But I don't know how many people want to do focal things and massive grids.  I think you can do star proxy things with roster as well. I have to look at it, because it's essentially rest of the same thing right so says here's the file name, and then I'm going to work with that and sort of return when I'm done.  So yeah that is that is relatively they work to sort of stars provides the from roster to roster conversions.  So you can do call call things there. It's convenient. In any case you don't have to sort of work out your band and time things.  Yeah, and the other thing is the conversion integration with SF, which is pretty natural because it's kind of that sense it starts extends SF for multi dimensional feature or feature time series, or multi dimensional feature data.  And then another questions.  So, yeah, so, as I said, there is this arm art down so I hope that some people are some people who actually managed to run the whole thing.  So in the meantime, yeah, good. Okay.  So in the meantime while I was talking it's sort of the run all button on our studio after installing a few bigger packages.  So when you don't manage then, you know, catch up with me in some breaks sometime.  I'll be here until Thursday morning.  I'm happy to discuss things are happy to also explore all the things you cannot do with with stars.
OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_0_11879.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_11880_11999.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_12000_12095.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_12096_13703.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_13704_14123.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_14124_15155.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_15156_16031.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_16032_18287.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_18288_19739.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_19740_24623.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_24624_26111.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_26112_26279.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_26280_26519.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_26520_27947.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_27948_29039.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_29040_29171.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_29172_29363.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_29364_29951.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_29952_30011.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_30012_30275.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_30276_31619.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_31620_31799.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_31800_31859.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_31860_32675.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_32676_33431.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_33432_33731.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_33732_33803.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_33804_34595.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_34596_34679.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_34680_36263.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_36264_36383.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_36384_39971.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_39972_40151.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_40152_41927.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_41928_43151.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_43152_43499.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_43500_43619.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_43620_44183.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_44184_44327.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_44328_44771.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_44772_44999.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_45000_45167.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_45168_47387.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_47388_48731.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_48732_49127.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_49128_53963.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_53964_54311.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_54312_55391.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_55392_55859.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_55860_56567.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_56568_78959.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_78960_82427.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_82428_83531.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_83532_83579.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_83580_84683.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_84684_84803.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_84804_87227.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_87228_91955.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_91956_95663.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_95664_100607.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_100608_102179.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_102180_103691.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_103692_103835.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_103836_104183.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_104184_105395.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_105396_107327.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_107328_108731.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_108732_113627.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_113628_113939.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_113940_114587.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_114588_120947.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_120948_122939.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_122940_124595.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_124596_126287.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_126288_128111.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_128112_131267.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_131268_133331.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_133332_140963.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_140964_141323.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_141324_141431.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_141432_141695.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_141696_143999.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_144000_145031.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_145032_145475.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_145476_148523.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_148524_148883.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_148884_149423.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_149424_151859.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_151860_153719.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_153720_155315.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_155316_155471.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_155472_155531.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_155532_156863.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_156864_166907.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_166908_167891.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_167892_167999.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_168000_168119.png OpenGeoHub_Summer_School_2020_Edzer_Pebesma-_Handling_and_analysing_vector_and_raster_data_cubes_wit_168120_171864.png