You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: _posts/rastertemplate/rastertemplate.Rmd
+37-23
Original file line number
Diff line number
Diff line change
@@ -1,7 +1,8 @@
1
1
---
2
2
title: "Raster Operations"
3
3
description: |
4
-
Learn how to obtain a template raster from the `{d6geodata}` package, expanding the raster in different ways, rasterizing vector data, and extracting information from a raster using vector data. The methods described involve the use of packages like `{terra}` and `{sf}`, and include examples and plots to illustrate the processes. Additionally, we show the use of the {tidyterra} package for plotting raster data with ggplot.
4
+
Learn how to 1) obtain a template raster from the `{d6geodata}` package, 2) expand the raster in different ways, 3) rasterize vector data, and 4) extract information from a raster using vector data.
5
+
The methods described involve the use of packages like `{terra}` and `{sf}`, and include examples and plots to illustrate the processes. Additionally, we show how to use the `{tidyterra}` package for plotting raster data with `{ggplot2}`.
5
6
6
7
preview: ./img/wiki/preview-rastertemplate.png
7
8
categories:
@@ -47,7 +48,7 @@ theme_set(
47
48
48
49
## Raster Template {#template}
49
50
50
-
First get the template raster. The origin is the extent of Berlin with the projection laea (EPSG:3035).
51
+
First get the template raster. The origin is the extent of Berlin with the projection LAEA (EPSG:3035).
51
52
52
53
```{r template}
53
54
temp <- d6geodata::temp_ras(1000) # resolution is 1000m
@@ -62,37 +63,41 @@ We load Berlin districts data set from our cloud to show the borders of Berlin.
62
63
berlin <- d6geodata::get_geodata(data_name = "districs_berlin_2022_poly_03035_gpkg",
63
64
path_to_cloud = "E:/PopDynCloud",
64
65
download_data = FALSE) %>%
65
-
st_union() %>% # using `st_union()` to combine
66
+
st_union() %>% # using `st_union()` to combine geometries
p1 <- d6geodata::plot_qualitative_map(temp) # you can use this build-in plot function from the `{d6geodata}` package
74
76
75
77
p1 + p_berlin
76
78
```
77
79
80
+
78
81
## Expand Raster {#expand}
79
82
80
83
### Expand On All Sides
81
84
82
-
If you want to expand the raster on all size (similar to a buffer around a polygon) you can just add cells to all sides. You can do this by adding an integer number in the y argument in the function `extend` of the `{terra}` package.
85
+
If you want to expand the raster (similar to a buffer around a polygon) you can just add cells to all sides. You can do this by adding an integer number in the y argument in the function `extend` of the `{terra}` package.
83
86
84
87
```{r}
88
+
# expand the map
85
89
temp_exp_all <- terra::extend(x = temp, # template raster as input
In this case we use a workaround by using a costum extent. We are adding on the left (western) side cells by using the xmin of the larger extend. Xmax, ymin, ymax are from the given extent.
112
+
In this case we use a workaround by using a custom extent. Here we are adding cells on the left (western) side by using the *xmin* of the larger extend. Xmax, ymin, ymax are from the given extent.
108
113
109
114
```{r}
110
115
temp_exp_1_sides <- terra::extend(x = temp,
@@ -123,19 +128,22 @@ p4 + p_berlin
123
128
Here you see a combination of all plots made to compare them directly.
plot_annotation(tag_levels = "1", tag_prefix = "(P", tag_suffix = ")") # add tag to every plot
129
135
```
130
136
131
137
## Rasterize Vector Data {#rasterize}
132
138
133
-
With `vect()` you can make a spatVector (a kind of shapfile) out of an extent object.
139
+
With `vect()` you can make a spatVector (similar to a shapefile) from an extent object.
134
140
Afterwards, we can use this spatVector to rasterize the data. This works similar with an sf object.
135
141
136
142
```{r plot}
137
-
temp_vec <- terra::vect(ext(temp), crs(temp)) # spatVector out of smalest raster extent
143
+
# create spatVector
144
+
temp_vec <- terra::vect(ext(temp), crs(temp)) # spatVector out of smallest raster extent
138
145
146
+
# rasterize spatVector
139
147
temp_rast <- rasterize(temp_vec, # vector data
140
148
temp_exp_all %>% disagg(10),
141
149
# raster data with larger extent, disaggregated to 100m
@@ -144,28 +152,34 @@ temp_rast <- rasterize(temp_vec, # vector data
144
152
temp_rast
145
153
```
146
154
155
+
147
156
### Rasterized Plot
148
157
158
+
View the results.
159
+
149
160
```{r}
150
161
d6geodata::plot_qualitative_map(temp_rast)
151
162
```
152
163
153
164
## Extract With Vector Data {#extract}
154
165
155
-
Imagine you have a track of your animal and you want to see, which type of information lays in this track. In this case, you can use your line and extract the values from a raster below. There are several ways to this but the easiest way is to use the `extract()` function form the `{terra}` package to get the values laying under the vector data.
166
+
Imagine you have a track of an animal and you want to know what information lays in this track. In this case, you can use your line and extract the values from a raster below. There are several ways to this but the easiest way is to use the `extract()` function form the `{terra}` package to get the values laying under the vector data.
156
167
157
168
For this, you can use any vector data (points, lines, polygons) with the same projection and extent to extract values from a raster. The `extract()` function uses the raster as first input and the vector data as second. Depending on the type of data you are using you can define a function for the extraction like min, max or mean. 'Cells' gives you the call ID, 'xy' the coordinates and 'ID' the row number of your extraction.
169
+
As mentioned, it is important that both the vector and the raster have the same projection.
158
170
159
171
### Line For Extraction
160
172
161
-
For this we just creating a dummy line with coordinates laying within the raster
173
+
For this we are just creating a dummy line with coordinates laying within the raster
162
174
163
175
```{r}
176
+
# create points
164
177
pts <- tibble(x = c(4528790, 4576890), # a point out of two coordinates from the raster
165
178
y = c(3271740, 3272140))
166
179
180
+
# create line from points
167
181
dummy_line <- st_as_sf(pts, # a line made out of the point data
168
-
coords = c("x", "y"),
182
+
coords = c("x", "y"), # specify columns of xy coordinates
169
183
crs = st_crs(temp)$wkt, # using the same projection as for the template raster
170
184
sf_column_name = "geometry",
171
185
remove = FALSE) %>%
@@ -177,8 +191,8 @@ dummy_line <- st_as_sf(pts, # a line made out of the point data
p5 <- d6geodata::plot_qualitative_map(temp_rast) + # first plot raster
195
+
geom_sf(data = dummy_line) # plot line on top of the raster
182
196
183
197
p5
184
198
@@ -193,31 +207,31 @@ ext_tab <- terra::extract(x = temp_rast, # template raster as first input
193
207
y = dummy_line, # dummy line as second input
194
208
fun=NULL, # here you can set a function to summarise the data directly
195
209
method="simple",
196
-
cells=FALSE, # here you can get the cell number as well
197
-
xy=FALSE, # here the xy coordinates
210
+
cells=FALSE, # here you can get the cell number as well if cell = TRUE
211
+
xy=FALSE, # here the xy coordinates if xy = TRUE
198
212
ID=TRUE # and the id in front of each row
199
213
)
200
214
```
201
215
202
-
As you already can see on the plot, the largest part of the line lays within the purple area (1). This function can be used in cases like you need to know in which landuse class your track lays.
216
+
As you already can see on the plot, the largest part of the line lays within the purple area (value = 1). This function can be used to know on which land use class your track lays, for example.
203
217
204
218
```{r}
205
219
table(ext_tab$layer)
206
220
```
207
221
208
-
## Visualize data with tidyterra {#tidyterra}
222
+
## Visualize Data With Tidyterra {#tidyterra}
209
223
210
-
There are several ways to plot raster data with ggplot and one relatively new way is to do it with the package `{tidyterra}`. It has some well developed `{ggplot}` functions, but can sometimes not be combined with the base functions of `{ggplot}`. If you use this pacakge you may have to stay within the package functions.
224
+
There are several ways to plot raster data with ggplot and one relatively new way is to do it with the package `{tidyterra}`. It has some well developed `{ggplot}` functions, but it cannot always be combined with the base functions of `{ggplot}`. If you use this package you may have to stay within the package functions.
211
225
212
226
213
227
```{r}
214
228
library(tidyterra)
215
229
216
230
ggplot() +
217
-
geom_spatraster(data = temp_rast) +
218
-
p_berlin
219
-
231
+
geom_spatraster(data = temp_rast) + # use geom_spatraster() to plot the raster data with {tidyterra}
0 commit comments