1

##Week 6

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

Warning: package 'sf' was built under R version 4.3.3

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; 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("tidyr")
library("terra")

Warning: package 'terra' was built under R version 4.3.3

terra 1.8.29

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract

library("stars")

Warning: package 'stars' was built under R version 4.3.3

Loading required package: abind

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.

##Task 1: Import and visualize spatial data

###Import #Feldaufnahmen

feldaufnahmen<- read_sf("Feldaufnahmen_Fanel.gpkg")

st_layers("Feldaufnahmen_Fanel.gpkg") #One Layer

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

feldaufnahmen

Simple feature collection with 975 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2568099 ymin: 1199766 xmax: 2578824 ymax: 1207836
Projected CRS: CH1903+ / LV95
# A tibble: 975 × 3
   FieldID Frucht                                                           geom
     <dbl> <chr>                                                   <POLYGON [m]>
 1       1 Roggen ((2570914 1202743, 2570917 1202749, 2570985 1202833, 2571294 …
 2       0 <NA>   ((2570893 1202758, 2570893 1202758, 2570959 1202845, 2570985 …
 3       0 <NA>   ((2570868 1202776, 2570872 1202781, 2570913 1202828, 2570946 …
 4       2 Wiese  ((2570882 1203234, 2570641 1202974, 2570630 1202983, 2570606 …
 5       3 Weide  ((2570249 1203116, 2570371 1203328, 2570481 1203197, 2570390 …
 6       5 Weide  ((2570378 1203320, 2570466 1203436, 2570552 1203289, 2570481 …
 7       6 Weide  ((2570466 1203436, 2570572 1203495, 2570659 1203433, 2570659 …
 8       4 Weide  ((2569706 1203278, 2569706 1203342, 2570199 1203198, 2570223 …
 9       7 Wiese  ((2570804 1203310, 2570805 1203312, 2570900 1203608, 2571208 …
10       0 Wald   ((2571004 1202990, 2571041 1203029, 2571073 1203003, 2571035 …
# ℹ 965 more rows

#simple plotting, shows already quiet a lot. 
st_crs(feldaufnahmen) #CH1903+ / LV95

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(feldaufnahmen["Frucht"])
plot(st_geometry(wildschwein_BE), add = TRUE)

1. What information does the dataset contain? The data set contrain information about field use, e.g. Roggen, Wiese, Weise.

2. What is the geometry type of the dataset (possible types are: Point, Lines and Polygons)? Polygons

3. What are the data types of the other columns? ID - Num and Frucht is a character (or ev. Factor)

4. What is the coordinate system of the dataset? Swiss Coordinate System CH1903+ / LV95

##Task 2: Annotate Trajectories from vector data

Join the two data set with st_join

joined <- st_join(wildschwein_BE,feldaufnahmen)

#st_joined is by default a left join. If you want a right join, you need to turn around the parameter (x,y)

##Task 3: Explore annotated trajectories

library("lubridate")

Attaching package: 'lubridate'

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

    intersect, union

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

    date, intersect, setdiff, union

library("forcats")
library("forcats")

wildschwein_smry<- joined |> 
  st_drop_geometry() |> 
  mutate(
    hour=hour(round_date(DatetimeUTC,"hour")),
    Frucht=fct_lump(Frucht,3)) |> 
  group_by(TierName,hour,Frucht) |> 
  summarise(n=n()) |> 
  group_by(TierName,hour) |> 
  mutate(perc = n/sum(n)*100)

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

wildschwein_smry |> 
  ggplot(aes(hour, perc, fill=Frucht))+
  geom_col()+
  facet_wrap(~TierName)+
  scale_y_continuous(
    labels = scales::label_number(suffix = "%"))+
  labs(
    x = "Time (rounded to the nearest hour)",
    y = "Percentage",
    title = "Percentages of samples in a given crop per hour",
    subtitle = "(Only showing the most common categories)"
  )

##Task 4: Import and visualize vegetationindex (raster data)

library("terra")

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

##Task 5 : Annotate Trajectories from raster data

wildboar_veg <- extract(veg, wildschwein_BE)

nrow(wildboar_veg)

[1] 51246

nrow(wildschwein_BE)

[1] 51246

wildschwein_BE$veg <- wildboar_veg$vegetationshoehe_LFI

wildschwein_BE

Simple feature collection with 51246 features and 7 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 2568153 ymin: 1202306 xmax: 2575154 ymax: 1207609
Projected CRS: CH1903+ / LV95
# A tibble: 51,246 × 8
   TierID TierName CollarID DatetimeUTC                E        N
 * <chr>  <chr>       <dbl> <dttm>                 <dbl>    <dbl>
 1 002A   Sabi        12275 2014-08-22 21:00:12 2570409. 1204752.
 2 002A   Sabi        12275 2014-08-22 21:15:16 2570402. 1204863.
 3 002A   Sabi        12275 2014-08-22 21:30:43 2570394. 1204826.
 4 002A   Sabi        12275 2014-08-22 21:46:07 2570379. 1204817.
 5 002A   Sabi        12275 2014-08-22 22:00:22 2570390. 1204818.
 6 002A   Sabi        12275 2014-08-22 22:15:10 2570390. 1204825.
 7 002A   Sabi        12275 2014-08-22 22:30:13 2570387. 1204831.
 8 002A   Sabi        12275 2014-08-22 22:45:11 2570381. 1204840.
 9 002A   Sabi        12275 2014-08-22 23:00:27 2570316. 1204935.
10 002A   Sabi        12275 2014-08-22 23:15:41 2570393. 1204815.
# ℹ 51,236 more rows
# ℹ 2 more variables: geometry <POINT [m]>, veg <dbl>

wildschwein_smry_veg <- wildschwein_BE |>
  st_drop_geometry() |>
  mutate(
    hour = hour(round_date(DatetimeUTC, "hour"))
  ) |>
  group_by(TierName, hour) |>
  summarise(
    veg_avg = mean(veg, na.rm = TRUE),
    .groups = "drop"
  )

wildschwein_smry_veg |>
  ggplot(aes(x = hour, y = veg_avg, fill = TierName)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~TierName) +
  labs(
    x = "Time",
    y = "Mean vegetation height [m]",
    title = "Vegetation height per animal and hour"
  )