Grassland continuity in sampled areas

Fist step is to read in coordinates of sampled sites, the dataset actually has 3 different datasets in it so it will have to be wrangled, you can download the excel sheet in the following link

coordinates <- read_excel("ExampleData/All Fields GPS coordinates.xlsx", skip = 1)

Coordinates_a <- coordinates[, 1:3] |>
    janitor::clean_names()
colnames(Coordinates_a) <- stringr::str_sub(colnames(Coordinates_a), end = -3)

Coordinates_a$type <- "Grassland"

Coordinates_b <- coordinates[, 5:7] |>
    janitor::clean_names()

colnames(Coordinates_b) <- stringr::str_sub(colnames(Coordinates_b), end = -3)

Coordinates_b$type <- "Conservation agriculture"

Coordinates_c <- coordinates[, 9:11] |>
    janitor::clean_names()

colnames(Coordinates_c)[1] <- stringr::str_sub(colnames(Coordinates_c)[1],
    end = -3)

colnames(Coordinates_c)[2:3] <- stringr::str_sub(colnames(Coordinates_c)[2:3],
    end = -4)

Coordinates_c$type <- "Conventional agriculture"

coordinates <- list(Coordinates_a, Coordinates_b, Coordinates_c) |>
    purrr::reduce(dplyr::bind_rows) |>
    tidyr::separate(coordinates_latitude_longitude, into = c("lat", "lon"),
        sep = ",") |>
    dplyr::mutate(lat = as.numeric(lat), lon = as.numeric(lon)) |>
    dplyr::mutate(location_code = stringr::str_replace_all(location_code,
        "Ø", "oe")) |>
    dplyr::mutate(location_code = stringr::str_replace_all(location_code,
        "Å", "aa")) |>
    dplyr::mutate(location_code = stringr::str_replace_all(location_code,
        "Æ", "ae")) |>
    dplyr::mutate(location_code = janitor::make_clean_names(location_code,
        allow_dupes = T)) |>
    dplyr::mutate(location_code = stringr::str_remove_all(location_code, "ca_"),
        location_code = stringr::str_remove_all(location_code, "k_"), location_code = stringr::str_remove_all(location_code,
            "_"))

coordinates <- coordinates[complete.cases(coordinates), ]

the final dataset can be downloaded here

Getting habitat proportion with given locations

First we will read in the spatRaster used for the habitat proportion calculation, the full resolution map can be downloaded in this link

data("Landuse_DK")
Landuse <- terra::unwrap(Landuse_DK)

Extract proportion for the points

coordinates_sf <- terra::vect(coordinates, geom = c("lon", "lat"), crs = "epsg:4326") |>
    terra::project(terra::crs(Landuse))

The landuse and points can be seen in the following figure

Extract proportion of grasland from points

GrasslandProp <- SpatioTemporalCont::summarise_polygons(Rast = Landuse, Polygons = coordinates_sf,
    Vars = "dry nature", dist = 2000, type = "Both")

dry nature
0.0312000
0.1876997
0.0043764
0.0000000
0.0175439
0.0276243
0.0306452
0.0079618
0.0192308
0.0087789
0.0091827
0.0520416
0.0116580
0.0335463
0.0240577
0.0176000
0.0076220
0.0098684
0.0127694
0.0333988
0.0069541
0.0179211
0.0240192
0.0971338
0.0122549
0.0207171
0.0713202
0.0192000
0.0341615
0.0511591
0.0000000
0.0087789
0.0059172
0.0287770
0.0064655
0.0120000
0.0422648
0.0311751
0.0255795
0.0664000
0.0486056
0.0447642
0.0641026
0.0535572
0.0047923
0.0638468
0.0215483
0.0307329
0.0087789
0.0239425
0.0120000
0.0112269
0.0160385
0.0275150
0.0645933
0.0368295
0.0247604
0.0000000
0.0023962
0.0024019
0.0000000
0.0000000
0.0000000
0.0020101
0.0000000
0.0016051
0.0000000
0.0000000
0.0000000
0.0111732
0.0281350
0.0056045
0.0055777
0.0088000
0.0023981
0.0000000
0.0079618
0.0000000
0.0962521
0.0056000
0.0000000
0.0084666
0.0128102
0.0008000
0.0016000

You can download the proportion table here

Importing genetic diversity

We will import the data diversity, you can download the raw file here

test_data_genetic_diversity_workshop2024 <- read_table("ExampleData/test_data_genetic_diversity_workshop2024.txt") |>
    dplyr::mutate(Site = janitor::make_clean_names(Site, allow_dupes = T),
        Site = stringr::str_remove_all(Site, "_"), Site = stringr::str_replace_all(Site,
            "lss", "les")) |>
    dplyr::select(Site, pi) |>
    dplyr::rename(location_code = Site) |>
    dplyr::full_join(coordinates)

Now we have joined most of the sites, and we will now filter to only sites with genetic diversity values, that are grasslands and have coordinates

clean_data <- test_data_genetic_diversity_workshop2024 |>
    dplyr::filter(!is.na(pi), type == "Grassland", !is.na(lat))

Lets look at genetic diversity in the map

clean_data_sf <- terra::vect(clean_data, geom = c("lon", "lat"), crs = "epsg:4326") |>
    terra::project(terra::crs(Landuse))

Genetic distance

First we generate a table of ids for the site

SiteID <- data.frame(location_code = c("aaRJ", "BIJ", "DoeJ", "DSJ", "FUR",
    "GEJ", "GUS", "HHJ", "JEJ", "JHJ", "KoeJ", "LVJ", "MSJ", "NOJ", "SBJ",
    "SKJ", "ULJ", "UTJ", "VAJ"), ID = 1:19) |>
    dplyr::mutate(location_code = stringr::str_to_lower(location_code)) |>
    dplyr::left_join(clean_data)

We also build a spatial version of this to get the spatial distance

SiteID_sf <- terra::vect(SiteID, geom = c("lon", "lat"), crs = "epsg:4326") |>
    terra::project(terra::crs(Landuse))

CoordSites <- terra::geom(SiteID_sf)[, c(3, 4)]

Now lets combine this with the actual Fst data

Fst <- read_table("ExampleData/EntNic_fst_data_old.tsv", col_names = FALSE) |>
    tidyr::separate(X1, into = c("from", "to")) |>
    dplyr::rename(Fst = X2) |>
    dplyr::mutate_all(as.numeric) |>
    dplyr::group_by(from, to) |>
    dplyr::summarise(Fst = mean(Fst))

And now add the euclidean distance

euc_dist <- distance(SiteID_sf) |>
    as.matrix() |>
    GeNetIt::dmatrix.df() |>
    right_join(Fst)

An now we can plot the distances among them

Building transition distances

We will build transistion distances based on 2500 meters and 500 meters for grasslands.

Grassland500 <- SpatioTemporalCont::calculate_prop(Rast = Landuse, Radius = 500,
    Vars = c("dry nature", "Forest", "wet nature", "Agriculture"))

Grassland500 <- Grassland500["dry nature"] * 2 + Grassland500["wet nature"] +
    Grassland500["Forest"] * 0.8 + Grassland500["Agriculture"] * 0.2

Grassland500 <- Grassland500 + 0.001
Grassland500[is.na(Grassland500)] <- 0.001

Grassland500 <- terra::aggregate(Grassland500, fact = 5, cores = 4)

Grassland2500 <- SpatioTemporalCont::calculate_prop(Rast = Landuse, Radius = 2500,
    Vars = c("dry nature", "Forest", "wet nature", "Agriculture"))

## 
|---------|---------|---------|---------|
=========================================

Grassland2500 <- Grassland2500["dry nature"] * 2 + Grassland2500["wet nature"] +
    Grassland2500["Forest"] * 0.8 + Grassland2500["Agriculture"] * 0.2

Grassland2500 <- Grassland2500 + 0.001
Grassland2500[is.na(Grassland2500)] <- 0.001

Grassland2500 <- terra::aggregate(Grassland2500, fact = 5, cores = 4)

We can plot both together

To build the conductivity layers we use gdistance

tr.cost500 <- gdistance::transition(raster::raster(Grassland500), transitionFunction = mean,
    directions = 4)
tr.cost500 <- gdistance::geoCorrection(tr.cost500, type = "c", multpl = FALSE)


saveRDS(tr.cost500, "tr.cost500.rds")

tr.cost2500 <- gdistance::transition(raster::raster(Grassland2500), transitionFunction = mean,
    directions = 4)
tr.cost2500 <- gdistance::geoCorrection(tr.cost2500, type = "c", multpl = FALSE)

saveRDS(tr.cost2500, "tr.cost2500.rds")

Check some fast graphical outputs first with euclidean distance

Neighbours <- terra::delaunay(clean_data_sf, as.lines = T)
terra::writeVector(Neighbours, "Neighbours.shp", overwrite = T)

and then with shortest paths

Lines500 <- list()

for (i in 1:nrow(Neighbours)) {
    Temp <- as.matrix(as.data.frame(terra::intersect(clean_data_sf, Neighbours[i,
        ]), geom = "XY")[c("x", "y")])
    Lines500[[i]] <- terra::vect(sf::st_as_sf(gdistance::shortestPath(tr.cost500,
        origin = Temp[1, ], goal = Temp[2, ], output = "SpatialLines")))
    terra::crs(Lines500[[i]]) <- terra::crs(Neighbours)
}

Lines2500 <- list()

for (i in 1:nrow(Neighbours)) {
    Temp <- as.matrix(as.data.frame(terra::intersect(clean_data_sf, Neighbours[i,
        ]), geom = "XY")[c("x", "y")])
    Lines2500[[i]] <- terra::vect(sf::st_as_sf(gdistance::shortestPath(tr.cost2500,
        origin = Temp[1, ], goal = Temp[2, ], output = "SpatialLines")))
    terra::crs(Lines2500[[i]]) <- terra::crs(Neighbours)
}

Lines500 <- Lines500 |>
    purrr::reduce(rbind)
terra::writeVector(Lines500, "Lines500.shp", overwrite = T)

Lines2500 <- Lines2500 |>
    purrr::reduce(rbind)
terra::writeVector(Lines2500, "Lines2500.shp", overwrite = T)

Now we plot the shortest path, vs euclidean

Landuse <- terra::unwrap(Landuse_DK)
Lines500 <- terra::vect("Lines500.shp")
Lines2500 <- terra::vect("Lines2500.shp")
Neighbours <- terra::vect("Neighbours.shp")


ggplot() + geom_spatraster(data = Landuse) + geom_spatvector(data = Neighbours) +
    geom_spatvector(data = Lines500, color = "red") + geom_spatvector(data = Lines2500,
    color = "blue")

Landuse <- terra::unwrap(Landuse_DK)
euc_dist$least_cost500 <- NA
euc_dist$least_cost2500 <- NA

for (i in 1:nrow(euc_dist)) {
    Temp <- terra::vect(sf::st_as_sf(gdistance::shortestPath(tr.cost500, origin = CoordSites[euc_dist$from[i],
        ], goal = CoordSites[euc_dist$to[i], ], output = "SpatialLines")))
    terra::crs(Temp) <- terra::crs(Landuse)
    euc_dist$least_cost500[i] <- terra::perim(Temp)

    Temp <- terra::vect(sf::st_as_sf(gdistance::shortestPath(tr.cost2500, origin = CoordSites[euc_dist$from[i],
        ], goal = CoordSites[euc_dist$to[i], ], output = "SpatialLines")))
    terra::crs(Temp) <- terra::crs(Landuse)
    euc_dist$least_cost2500[i] <- terra::perim(Temp)
}

all_dist <- euc_dist

Now we can show the plot of the different distances

Buidling the mesh

First we will take a look at a first option of a mesh

Meshcoords <- terra::geom(clean_data_sf) |>
    as.data.frame()

bnd <- INLA::inla.nonconvex.hull(cbind(Meshcoords$x, Meshcoords$y), convex = -0.1)
mesh_inla <- fmesher::fm_mesh_2d(boundary = bnd, max.edge = c(8000, 50000))
mesh <- sdmTMB::make_mesh(Meshcoords, c("x", "y"), mesh = mesh_inla)

We should add a barrier, since we have the sea, for that we will download a polygon of Denmark, reproject it and transform it in to an sf object

DK <- geodata::gadm("Denmark", level = 0, path = getwd()) |>
    terra::project(Landuse) |>
    sf::st_as_sf()

And now we add it to the mesh

bspde <- sdmTMBextra::add_barrier_mesh(mesh, DK, range_fraction = 0.1, plot = FALSE)

Afternoon work

Derek Corcoran

2024-02-24

Packages used