Exercise_6

Task 1

library("readr")
library("sf")

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

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.

feld <- read_sf("Feldaufnahmen_Fanel.gpkg")


str(feld)

sf [975 × 3] (S3: sf/tbl_df/tbl/data.frame)
 $ FieldID: num [1:975] 1 0 0 2 3 5 6 4 7 0 ...
 $ Frucht : chr [1:975] "Roggen" NA NA "Wiese" ...
 $ geom   :sfc_POLYGON of length 975; first list element: List of 1
  ..$ : num [1:12, 1:2] 2570914 2570917 2570985 2571294 2571719 ...
  ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
 - attr(*, "sf_column")= chr "geom"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  ..- attr(*, "names")= chr [1:2] "FieldID" "Frucht"

st_crs(feld)

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 of initial line",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 on initial line",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]]

plot(st_geometry(feld))
plot(st_geometry(wildschwein_BE), add = TRUE, reset = FALSE)

What information does the dataset contain?
  ID of the field, type of fruit, polygon of field
What is the geometry type of the dataset (possible types are: Point, Lines and Polygons)?
  Polygon
What are the data types of the other columns?
  ID = num, Frucht = chr, geom = polygon
What is the coordinate system of the dataset?
  CH1903+ / LV95, 2056
  
  

Task 2

feldschwein <- st_join(wildschwein_BE, feld, join = st_within)

Task 3

library(ggplot2)
library(lubridate)

Attaching package: 'lubridate'

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

    date, intersect, setdiff, union

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

feldschwein2 <- feldschwein |> 
  st_drop_geometry() |> 
  mutate(
    hour = hour(round_date(DatetimeUTC, "hour"))
  ) |> 
  #count(Tiername, hour, Frucht)
  group_by(TierName, hour, Frucht) |> 
  summarise(n = n()) |> 
  group_by(TierName, hour) |> 
  mutate(perc = n/sum(n))

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

feldschwein2 |> 
  ggplot(aes(hour, perc, fill = Frucht)) +
  geom_col() +
  facet_wrap(~TierName)

Task 4

library(terra)

terra 1.7.78

library(tmap)

veg <- rast("vegetationshoehe_LFI.tif")

plot(veg)

tm_shape(veg) +
  tm_raster(style = "cont", palette = "viridis", title = "Heatmap") +
  tm_layout(legend.outside = TRUE)

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
  'tm_scale_continuous(<HERE>)'
[v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
SpatRaster object downsampled to 2753 by 3634 cells.

Task 5

boarveg <- extract(veg, wildschwein_BE) #y needs to be a SpatVector object

Warning: [extract] transforming vector data to the CRS of the raster

nrow(boarveg)

[1] 51246

nrow(wildschwein_BE)

[1] 51246

# Same amount of rows

wildschwein_BE$veg <- boarveg$vegetationshoehe_LFI

wildschwein_BE$veg_group <- cut(wildschwein_BE$veg, breaks = 5, labels = c("Very low", "Low", "Middle", "high", "Very high"), include.lowest = TRUE )

wildschwein_BE2 <- wildschwein_BE |> 
  st_drop_geometry() |> 
  mutate(
    hour = hour(round_date(DatetimeUTC, "hour"))
  ) |> 
  #count(Tiername, hour, veg_group)
  group_by(TierName, hour, veg_group) |> 
  summarise(n = n()) |> 
  group_by(TierName, hour) |> 
  mutate(perc = n/sum(n))

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

wildschwein_BE2 |> 
  ggplot(aes(hour, perc, fill = veg_group)) +
  geom_col() +
  facet_wrap(~TierName)