Up to this point I used reference levels (RL) drawn from my own personal and subjective intuition. There might be more objective ways to do this, or at least more precise and defendable ways. What we have decided on is to create a set of calibration maps, with varying amounts of trenching, and then several experts will look at these and judge individually how severe they believe a given density of trenches will be for the mire ecosystem. We need (minimum) three RLs. We define our RL based on the 7GR-GI variable in Halvorsen and Bratli (2019).
There is no way for us to know the ages of the trenches in our data. Therefore it is impossible to separate the realised from the future effects from trenching. We therefore combine these two aspects, and describe the ecosystem effect of trenching based on the future equilibrium state that will result as a consequence of the trenching that has happened. This approach is not ideal, and the ecosystem condition should preferably be describes as a snap shot of the condition as it is today. However, we see know other solution to this issue.
X0 is the variable value (i.e. the frequency of 10x10 grid cells with trenching inside each 100x100 cell) that corresponds to a completely destroyed ecosystem. A further deterioration of the variable value should now result in any further deterioration of ecosystem condition. We have defined this as a 7GR-GI value of 5, meaning that the ecosystem is set on a succession towards a non-wetland type, usually forest.
X100 is the RL at the reference condition (RC; pristine condition in our case). This RL is quite intuitive and is set to zero (7GR-GI = 1). Note that the trenches data has been filtered to remove very small polygons that are more likely to be artefacts than real trenches.
X60 represents the variable value at the boundary between good and deteriorated ecosystem condition. We have decided to anchor this RL theoretically to the point where trenching is so severe that the species composition in the mire ecosystem has changes so much as to correspond to a one-step change along one of the dominant environmental gradient (LKMs; Halvorsen and Bratli (2019)). This equates to 7GR-GI = 3.
In practice we map 7GR-GI = 1 to a value between 0.8 and 1 on the indicator scale, 7GR-GI = 2 to values between 0.6 and 0.4, 7GR-GI = 3 to values between 0.4 and 0.6, 7GR-GI = 4 to values between 0.2 and 0.4, and 7GR-GI = 5 to values between 0.0 and 0.2,
7GR-GI is meant to be describes for whole hydrological units, but here we instead use it for individual grid cells. This means that a grid cell can be described as having no trenching, and this a favorable indicator value, even though the mire as a whole is severely affected by trenching arising from neighboring grid cells. This implies that our indicator could end up to be to conservative. However, we think the size of our grid cells (100x100m) is a relevant scale for describing trenching effects, that are known to often extend tens of meters from a large trench, but rarely hundreds of meters.
To set the RLs for our variable, we could potentially use empirical data, such as nature type mapping data where mire polygons are assign 7GR-GI values in the field. However, this field protocol does not include mapping of more common mire types, and although our study areas is mapped using this protocol, no mire polygons have been delineated.
6.1 Structured expert elicitation
Rather than using empirical data, we decided to use a kind of structured expert elicitation. Our expert panel consisted of four ecologist with extensive and relevant field experience. First we agreed on the theoretical definition of the RLs (se above). Then we looked at a set of 100x100m grid cells with varying densities of trenches (@fig-calplots). We then scored these cases individually, but following a common discussion, giving each case a precise indicator value. Our interpretation of the indicator scale was aided by both the five condition classes in the water framework directive (high, good, moderate, low, and poor), as well as the five steps on the 7GR-GI metric. This approach mean that we had the same understanding of the indicator scale and the general idea behind the exercise, but it also meant that our answers were influenced by the opinions of the other experts.
Set up Google API, download and save satelite images
# I used this to save the API key to my syst env# ggmap::register_google(key = "", write = TRUE)# example of getting sat image for the entire area# get center coordinatee<-st_as_sfc(st_bbox(grid100_freq4))|>st_transform(4326)|>st_centroid()|>st_coordinates()|>as.numeric()names(e)<-c("lon", "lat")satImg<-ggmap::get_googlemap( center =e, zoom =14, maptype ="satellite")#ggmap(satImg)#saveRDS(satImg, "data/satImg.rds")
Get some grid cells with varying amount of trenching
# Transform to 4326 to match ggmap objectscells<-grid100_freq4|># requre at least half a grid cellfilter(area>set_units(5000, m2))|>mutate(tens =round(freq, -1))|>group_by(tens)|>slice_max(freq, n =2, with_ties =F)|>ungroup()|>st_transform(4326)cells_c<-cells|>st_centroid()|>st_coordinates()|>as_tibble()|>rename("lon"="X", "lat"="Y")cells<-cells|>bind_cols(cells_c)
Function for downloading static google map
myGgmap<-function(x, zoom=18){ggmap::get_googlemap( center =c(lon =x$lon, lat =x$lat), zoom =zoom, maptype ="satellite")}# Hackish function to fix bbox in the ggmap object# https://stackoverflow.com/questions/47749078/how-to-put-a-geom-sf-produced-map-on-top-of-a-ggmap-produced-rasterggmap_bbox<-function(map){if(!inherits(map, "ggmap"))stop("map must be a ggmap object")# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, # and set the names to what sf::st_bbox expects:map_bbox<-setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))# Coonvert the bbox to an sf polygon, transform it to 3857, # and convert back to a bbox (convoluted, but it works)bbox_3857<-st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs =4326)), 3857))# Overwrite the bbox of the ggmap object with the transformed coordinates attr(map, "bb")$ll.lat<-bbox_3857["ymin"]attr(map, "bb")$ll.lon<-bbox_3857["xmin"]attr(map, "bb")$ur.lat<-bbox_3857["ymax"]attr(map, "bb")$ur.lon<-bbox_3857["xmax"]map}
# convert grid cell df to 3857cells2<-st_transform(cells, 3857)lev_crop2<-st_transform(lev_crop, 3857)for(iin1:nrow(cells2)){dat<-cells2[i, ]suppressMessages(without<-ggmap(ggmap =satimg[[i]])+coord_sf(crs =st_crs(3857))+geom_sf( data =dat, inherit.aes =F, fill =NA, color ="yellow", size =2))suppressMessages(with<-ggmap(ggmap =satimg[[i]])+coord_sf(crs =st_crs(3857))+geom_sf(data =dat, inherit.aes =F, fill =NA, color ="yellow", size=2)+geom_sf(data =lev_crop, inherit.aes =F, fill ="orange", color ="orange")+ggtitle(paste("Trenching frequency = ", round(dat$freq, 1),"%")))assign( x =paste0("ID", dat$id),ggarrange(with, without))}
We scored indicator values for the pictures above, and found that a lot of change was happening in the range between 6 and 17%. We therefore made these extra figures below to get some more data points for this range of trenching frequencies. Since we now had a good internal harmonized understanding of the task, we decided to score these cases indivudually, not knowing the answers or reasoning from the other experts.
Get some grid cells with varying amount of trenching between 6 and 17%
# Transform to 4326 to match ggmap objectscells_fine<-grid100_freq4|># requre at least half a grid cellfilter(area>set_units(5000, m2),between(freq, 6, 17))|>st_transform(4326)cells_c_fine<-cells_fine|>st_centroid()|>st_coordinates()|>as_tibble()|>rename("lon"="X", "lat"="Y")cells_fine<-cells_fine|>bind_cols(cells_c_fine)|>filter(!id%in%cells$id)
# convert grid cell df to 3857cells_fine2<-st_transform(cells_fine, 3857)#lev_crop2 <- st_transform(lev_crop, 3857)for(iin1:nrow(cells_fine)){dat<-cells_fine2[i, ]suppressMessages(without<-ggmap(ggmap =satimg2[[i]])+coord_sf(crs =st_crs(3857))+geom_sf( data =dat, inherit.aes =F, fill =NA, color ="yellow", size =2))suppressMessages(with<-ggmap(ggmap =satimg2[[i]])+coord_sf(crs =st_crs(3857))+geom_sf(data =dat, inherit.aes =F, fill =NA, color ="yellow", size=2)+geom_sf(data =lev_crop, inherit.aes =F, fill ="orange", color ="orange")+ggtitle(paste("Trenching frequency = ", round(dat$freq, 1),"%")))assign( x =paste0("ID", dat$id),ggarrange(with, without))}
Halvorsen, R., and Bratli, H. 2019. Dokumentasjon av NiN versjon 2.2 tilrettelagt for praktisk naturkartlegging: Utvalgte variabler fra beskrivelsessystemet. natur i norge, artikkel 11 (versjon 2.2.0).
