parameters_valugaps.R

library(raster)
library(terra)
library(sf)
library("rnaturalearth")
library("rnaturalearthdata")
library("exactextractr")


designpath<- "Projects/ValuGaps/Designs/"





hnv <- rast("Projects/ValuGaps/data/gis/hnv_germany.tif")

# Load in HNV raster data (preprocessed)


# Load BFN protected areas data

protected_areas <- readRDS("Projects/ValuGaps/data/gis/protected_bfn_raster.rds")

##############################################################################################
# Load in Protected areas data (preprocessed) (protected planet data, comment out if use bfn data)

# protected_areas_all <- readRDS("Projects/ValuGaps/data/gis/germany_protected.rds")
# levels(protected_areas_all)
# 
# categories_to_drop <- c(1,3, 4, 6, 8, 9, 11)   # define categories that should not be included in protected areas
# # keeps nationalparks, nature conservation areas and UNESCO biosphere reserve
# 
# protected_areas <- protected_areas_all
# protected_areas[protected_areas %in% categories_to_drop] <-NA    # Make binary categories
# 
# names(protected_areas) <- "value"
##############################################################################################


# Code protected areas as zero/one 
protected_areas <- ifel(is.na(protected_areas), 0, 1)



# Load in shapefile of germany
germany_sf <- read_sf("Projects/ValuGaps/data/gis/WDPA_WDOECM_Feb2023_Public_DEU_shp/gadm41_DEU_shp", "gadm41_DEU_0") # no county borders
germany_sf_county <- read_sf("Projects/ValuGaps/data/gis/WDPA_WDOECM_Feb2023_Public_DEU_shp/gadm41_DEU_shp", "gadm41_DEU_2")

germany_geo <- st_transform(germany_sf$geometry, crs=crs(protected_areas)) # transform to raster crs
germany_geo_county <- st_transform(germany_sf_county$geometry, crs=crs(protected_areas))


#### Read in CORINE Land cover data ####

corine <- rast("Projects/ValuGaps/data/gis/corine_raster_2018/DATA/U2018_CLC2018_V2020_20u1.tif")
germany_geo_c <- st_transform(germany_sf_county$geometry, crs=crs(corine)) # used to crop corine data to Germany
corine <- crop(corine, germany_geo_c)
corine <-  mask(corine, vect(germany_geo_c))
levels(corine)

hnv_categories <- c(12, 14:22, 26:29, 32, 35:37) # see p.37 HNV EUROPE DATA Documentation
corine_hnv <- corine
corine_hnv[!(corine_hnv %in% hnv_categories)] <-NA    # Make binary categories
corine_hnv <- ifel(is.na(corine_hnv), 0, 1)
corine_hnv <- mask(corine_hnv, vect(germany_geo_c))
corine_hnv <- mask(corine_hnv, hnv)


#### Protected areas CORINE
protected_area_categories <- c(10, 12:41) # remove artificial areas except UGS and remove marine waters
corine_pa <- corine
corine_pa[!(corine_pa %in% protected_area_categories)] <-NA    # Make binary categories
corine_pa <- ifel(is.na(corine_pa), 0, 1)
corine_pa <- mask(corine_pa, vect(germany_geo_c))

protected_areas <- mask(protected_areas, vect(germany_geo_c))
####



#### Define simulation parameters ####

resps = 96  # number of respondents 
nosim=1  # number of simulations to run (about 500 is minimum)
#buffer_dist <- 15000

buffer_dist_ls <- list(10000, 15000, 25000, 50000) # define radius for the respondents in meter


# betacoefficients should not include "-"
basc = 0.5
bb = 1.2
bb2 = -0.016
bc= 1.4
bc2 = -0.0175
bp=-0.03

mani1 = expr(left_join(final_set, by="ID"))

manipulations = list(
#expr(across(c("HNV_on_corine", "Protected_on_corine"), ~replace_na(.,0)) ),
alt1.protected = expr(alt1.b+Protected_on_corine*100),
alt2.protected = expr(alt2.b+Protected_on_corine*100),
alt3.protected = expr(Protected_on_corine*100),
alt1.HNV       = expr(alt1.b+HNV_on_corine*100),
alt2.HNV       = expr(alt2.b+HNV_on_corine*100),
alt3.HNV       = expr(HNV_on_corine*100),
expr(across(.cols=matches("HNV|protected"),.fns =  ~.x^2, .names = "{.col}{'sq'}" ) )


)

#place your utility functions here

u<-list( u1= list(
  v1 =V.1~ basc + bb * alt1.protected +  bb2 * alt1.protectedsq + bc * alt1.HNV + bc2 * alt1.HNVsq + bp * alt1.p , 
  v2 =V.2~ basc + bb * alt2.protected +  bb2 * alt2.protectedsq + bc * alt2.HNV + bc2 * alt2.HNVsq + bp * alt2.p ,
  v3 =V.3~        bb * alt3.protected +  bb2 * alt3.protectedsq + bc * alt3.HNV + bc2 * alt3.HNVsq )
)



## this function will be called before the data is simulated
sou_gis <- function() {
  
  source("Projects/ValuGaps/code/GIS_data.R")
  
}