library(tidyverse)
library(sf)
library(terra) # to work with raster data
wildschwein_BE <- read_delim("wildschwein_BE_2056.csv", ",") |>
st_as_sf(coords = c("E", "N"), crs = 2056, remove = FALSE)
wildschwein_BE$TierID <- as.factor(wildschwein_BE$TierID)
wildschwein_BE$TierName <- as.factor (wildschwein_BE$TierName)
Feldaufnahmen <- read_sf("Feldaufnahmen_Fanel.gpkg")
Vegetationshoehe <- rast("vegetationshoehe_LFI.tif") # import RasterExercise 6
Load libraries and data into R
Tasks 1: Import and visualize spatial data
Since Feldaufnahmen_Fanel.gpkg is a vector dataset, you can import it using read_sf(). Explore this dataset in R to answer the following questions:
What information does the dataset contain?
Field ID, Frucht, geom
What is the geometry type of the dataset (possible types are: Point, Lines and Polygons)?
It is a polygon. CRS: CH1903+ / LV95
What are the data types of the other columns?
- FieldID: numeric
- Frucht: chr
# change Frucht into a Factor
Feldaufnahmen$Frucht <- as.factor(Feldaufnahmen$Frucht)
unique(Feldaufnahmen$Frucht) [1] Roggen <NA> Wiese Weide Wald
[6] Gerste Sonnenblumen Mais Chinaschilf Acker
[11] Weizen Kartoffeln Erbsen Lupinen Flugplatz
[16] Rueben Raps Obstplantage Zwiebeln Zucchetti
[21] Blumenbeet Kohl Karotten Spinat Bohnen
[26] Hafer Salat Lauch Kuerbis Gewaechshaus
[31] Gemuese Kohlrabi Spargel Fenchel Mangold
[36] Brache Sellerie Rhabarber Feuchtgebiet Brokkoli
[41] Getreide Baumschule Buntbrache Flachs Ackerbohnen
44 Levels: Acker Ackerbohnen Baumschule Blumenbeet Bohnen Brache ... ZwiebelnWhat is the coordinate system of the dataset?
st_crs(Feldaufnahmen)Coordinate Reference System:
User input: CH1903+ / LV95
wkt:
PROJCRS["CH1903+ / LV95",
BASEGEOGCRS["CH1903+",
DATUM["CH1903+",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4150]],
CONVERSION["Swiss Oblique Mercator 1995",
METHOD["Hotine Oblique Mercator (variant B)",
ID["EPSG",9815]],
PARAMETER["Latitude of projection centre",46.9524055555556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8811]],
PARAMETER["Longitude of projection centre",7.43958333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8812]],
PARAMETER["Azimuth at projection centre",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8813]],
PARAMETER["Angle from Rectified to Skew Grid",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8814]],
PARAMETER["Scale factor at projection centre",1,
SCALEUNIT["unity",1],
ID["EPSG",8815]],
PARAMETER["Easting at projection centre",2600000,
LENGTHUNIT["metre",1],
ID["EPSG",8816]],
PARAMETER["Northing at projection centre",1200000,
LENGTHUNIT["metre",1],
ID["EPSG",8817]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
AREA["Liechtenstein; Switzerland."],
BBOX[45.82,5.96,47.81,10.49]],
ID["EPSG",2056]]Task 2: Annotate Trajectories from vector data
We would like to know what crop was visited by which wild boar, and at what time. Since the crop data is most relevant in summer, filter your wild boar data to the months may to june first and save the output to a new variable.
# create new dataframe for wildboar movement only in May to June
min(wildschwein_BE$DatetimeUTC)[1] "2014-08-22 21:00:12 UTC"max(wildschwein_BE$DatetimeUTC)[1] "2015-07-27 11:00:14 UTC"from <- as.POSIXct("2015-05-01", tz = "UTC")
to <- as.POSIXct("2015-06-30", tz = "UTC")
wildschwein_MJ <- wildschwein_BE |>
filter(DatetimeUTC >= from,
DatetimeUTC <= to)Overlay the filtered dataset with your fanel data to verify the spatial overlap.
# use spatial join to overlp datasets
wildschwein_join <- st_join(wildschwein_MJ, Feldaufnahmen)To sematically annotate each wild boar location with crop information, you can use a spatial join with the function st_join(). Do this and explore your annotated dataset.
Task 3: Explore annotated trajectories
Think of ways you could visually explore the spatio-temporal patterns of wild boar in relation to the crops. In our example below we visualize the percentage of samples in a given crop per hour.
wildschwein_join <- wildschwein_join |>
mutate(time_rounded = round_date(DatetimeUTC, "hour")) |>
mutate(hour = hour(time_rounded))
wildschwein_join$Frucht <- as.character(wildschwein_join$Frucht)
# rows to keep unchanged
keep <- c("Rueben", "Gerste", "Feuchtgebiet", "Wald")
wildschwein_join$Frucht[!wildschwein_join$Frucht %in% keep] <- "other"
# change back to a Factor
wildschwein_join$Frucht <- as.factor(wildschwein_join$Frucht)
wildschwein_join %>%
group_by(TierName, hour, Frucht) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(TierName, hour) %>%
mutate(percent = n / sum(n) * 100) %>%
ggplot(aes(hour, percent, fill = Frucht)) +
geom_col(width = 1) +
labs(
x = "Hour of day",
y = "Percentage of samples",
fill = "Crop"
) +
theme_minimal() +
facet_wrap(~TierName)
Task 4: Import and visualize vegetationindex (raster data)
In terms of raster data, we have prepared the Vegetation Height Model provided by the Swiss National Forest Inventory (NFI). This dataset contains high resolution information (1x1 Meter) on the vegetation height, which is determined from the difference between the digital surface models DSM and the digital terrain model by swisstopo (swissAlti3D). Buildings are eliminated using a combination of the ground areas of the swisstopo topographic landscape model (TLM) and spectral information from the stereo aerial photos.
Import this dataset using the function rast from the package terra. Visualize the raster data using the base plot function and tmap. We do not recommend using ggplot2 in this case, as is very slow with raster data.
library(tmap)
plot(Vegetationshoehe)
tm_shape(Vegetationshoehe) +
tm_raster()
Task 5: Annotate Trajectories from raster data
Semantically annotate your wild boar locations with the vegetation index (similar as you did with the crop data in Task 2). Since you are annotating a vector dataset with information from a raster dataset, you cannot use st_join but need the function extract from the terra package. Read the help on the extract function to see what the function expects.
wildschwein_BE %>%
mutate(vegetationshoehe_LFI = terra::extract(Vegetationshoehe,
terra::vect(.)
)[,2]) %>%
select(-CollarID, -E, -N, - TierName)Simple feature collection with 51246 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2568153 ymin: 1202306 xmax: 2575154 ymax: 1207609
Projected CRS: CH1903+ / LV95
# A tibble: 51,246 × 4
TierID DatetimeUTC geometry vegetationshoehe_LFI
<fct> <dttm> <POINT [m]> <dbl>
1 002A 2014-08-22 21:00:12 (2570409 1204752) 20.1
2 002A 2014-08-22 21:15:16 (2570402 1204863) 23.9
3 002A 2014-08-22 21:30:43 (2570394 1204826) 25.0
4 002A 2014-08-22 21:46:07 (2570379 1204817) 21.6
5 002A 2014-08-22 22:00:22 (2570390 1204818) 15.7
6 002A 2014-08-22 22:15:10 (2570390 1204825) 23.8
7 002A 2014-08-22 22:30:13 (2570387 1204831) 25.1
8 002A 2014-08-22 22:45:11 (2570381 1204840) 24.9
9 002A 2014-08-22 23:00:27 (2570316 1204935) 29.9
10 002A 2014-08-22 23:15:41 (2570393 1204815) 21.5
# ℹ 51,236 more rowswildschwein_extract <- extract(Vegetationshoehe, wildschwein_BE, bind = TRUE) %>% # creates a spatVector
sf::st_as_sf() |> # converts the spatVector into a sf object again
dplyr::select(
TierID,
DatetimeUTC,
vegetationshoehe_LFI,
geometry
)You can now explore the spatio-temporal patterns of this new data, similarly to as we did in Figure 24.1.
wildschwein_extract <- wildschwein_extract |>
mutate(time_rounded = round_date(DatetimeUTC, "hour")) |>
mutate(hour = hour(time_rounded))
wildschwein_extract %>%
group_by(TierID, hour, vegetationshoehe_LFI) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(TierID, hour) %>%
mutate(percent = n / sum(n) * 100) %>%
ggplot(aes(hour, percent, fill = vegetationshoehe_LFI)) +
geom_col(width = 1) +
labs(
x = "Hour of day (rounded)",
y = "Percentage",
fill = "Vegetationshoehe"
) +
theme_minimal() +
facet_wrap(~TierID)