Week6

Week 6

Exercises

Libraries

library(readr)
library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(ggplot2)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(lubridate)

Attaching package: 'lubridate'

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

library(terra)

terra 1.8.29

library(tmap)
library(forcats)

Read file

wildschwein_BE <- read_delim("wildschwein_BE_2056.csv", ",") |>
    st_as_sf(coords = c("E", "N"), crs = 2056, remove = FALSE)

Rows: 51246 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): TierID, TierName
dbl  (3): CollarID, E, N
dttm (1): DatetimeUTC

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Tasks 1: Import and visualize spatial data

st_layers("Feldaufnahmen_Fanel.gpkg")

Driver: GPKG 
Available layers:
           layer_name geometry_type features fields       crs_name
1 Feldaufnahmen_Fanel       Polygon      975      2 CH1903+ / LV95

field <- read_sf("Feldaufnahmen_Fanel.gpkg")
st_crs(field)

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]]

Explore this dataset in R to answer the following questions:

What information does the dataset contain? One layer with infos about Feldaufnahmen_Fanel.

What is the geometry type of the dataset (possible types are: Point, Lines and Polygons)? it is a Polygon-Type, with 975 features

What are the data types of the other columns? its a character with Frucht as a string. Field ID is nummeric.

What is the coordinate system of the dataset? CH1903+ / LV95

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. -> wildschwein_summer

Overlay the filtered dataset with your fanel data to verify the spatial overlap.

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.

# filter for may to june
wildschwein_summer <- wildschwein_BE |> 
  filter(month(DatetimeUTC) >= 5 & month(DatetimeUTC) <= 6)

# Overlay and plot fanel data to wildboar
plot(field["Frucht"], reset = FALSE)
plot(st_geometry(wildschwein_summer), add = T)

# join
wildschwein_field <- st_join(wildschwein_summer, field)
plot(wildschwein_field)

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.

summary <- wildschwein_field |>
  st_drop_geometry() |>
  mutate(
    hour = hour(round_date(DatetimeUTC, "hour")),
    Frucht = fct_lump(Frucht, 5)) |>
  group_by(TierName, hour, Frucht) |>
  summarise(n = n()) |>
  group_by(TierName, hour) |>
  mutate(percent = n/sum(n))

`summarise()` has grouped output by 'TierName', 'hour'. You can override using
the `.groups` argument.

ggplot(summary, aes(hour, percent, fill = Frucht)) +
  geom_col() +
  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.

vegetation <- rast("vegetationshoehe_LFI.tif")
plot(vegetation)

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.

# extract
wildschwein_vegetation <- extract(vegetation, wildschwein_summer)

# add vegetationhoehe to wildschwein_summer
wildschwein_summer$vegetation <- wildschwein_vegetation$vegetationshoehe_LFI
hist(wildschwein_summer$vegetation)

wildschwein_summer

Simple feature collection with 15559 features and 7 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 2568174 ymin: 1203130 xmax: 2574285 ymax: 1207596
Projected CRS: CH1903+ / LV95
# A tibble: 15,559 × 8
   TierID TierName CollarID DatetimeUTC                E        N
 * <chr>  <chr>       <dbl> <dttm>                 <dbl>    <dbl>
 1 002A   Sabi        12275 2015-05-01 00:00:17 2570093. 1205256.
 2 002A   Sabi        12275 2015-05-01 00:15:25 2570093. 1205249.
 3 002A   Sabi        12275 2015-05-01 00:30:15 2570091. 1205253.
 4 002A   Sabi        12275 2015-05-01 00:45:15 2570059. 1205242.
 5 002A   Sabi        12275 2015-05-01 01:00:30 2570078. 1205246.
 6 002A   Sabi        12275 2015-05-01 01:15:43 2570096. 1205256.
 7 002A   Sabi        12275 2015-05-01 01:30:17 2570079. 1205247.
 8 002A   Sabi        12275 2015-05-01 01:45:18 2570094. 1205248.
 9 002A   Sabi        12275 2015-05-01 02:00:45 2570044. 1205212.
10 002A   Sabi        12275 2015-05-01 02:15:22 2570025. 1205217.
# ℹ 15,549 more rows
# ℹ 2 more variables: geometry <POINT [m]>, vegetation <dbl>

veg_avg <- wildschwein_summer |>
  st_drop_geometry() |>
  group_by(TierName) |>
  summarise(avg_lfi = mean(vegetation, na.rm = TRUE))

buckets <- c(0, 1, 10, 20, 30, 40)
bucketed = cut(wildschwein_summer$vegetation, breaks = buckets)
bucketed_all <- ifelse(is.na(bucketed), "N/A", as.character(bucketed))
wildschwein_summer$veg_bucket <- bucketed_all

veg_monthly_lfi <- wildschwein_summer |>
  st_drop_geometry() |>
  mutate(month = month(DatetimeUTC)) |>
  group_by(TierName, month, veg_bucket) |>
  summarise(n = n()) |>
  group_by(TierName, month) |>
  mutate(percent = n/sum(n))

`summarise()` has grouped output by 'TierName', 'month'. You can override using
the `.groups` argument.

ggplot(data = veg_monthly_lfi, aes(x = as.factor(month), y = percent, fill = veg_bucket)) +
  geom_col() + 
  facet_wrap(~TierName) +
   labs(x = "Month", y = "Distribution of vegetation height")