working version with small changes

Projects/ValuGaps/code/GIS_data.R

-# Functions to compute shares 
-
-sum_cover <- function(x){
-  list(x %>%
-         group_by(value) %>%
-         summarize(total_area = sum(coverage_area, na.rm=T)) %>%
-         mutate(proportion = total_area/sum(total_area[1:2])))
-  
-}
-
-process_data <- function(data, name_prefix, cover_data) {
-  result <- exact_extract(data, buffer_pts, fun=sum_cover, coverage_area=T, summarize_df=T)
-  col_name <- paste(name_prefix)
-  result <- as.data.frame(sapply(result, function(x) ifelse(length(x$proportion) >= 2, x$proportion[2], 0)))
-  colnames(result) <- col_name
-  result <- result %>% mutate(ID = 1:resps) %>% dplyr::select(ID, everything())
-  return(result)
-}
-
-
-##### Load in GIS data ####
-
-## Draw random points without incorporating actual pop density (old)
-# random_points <- st_sample(germany_geo_county, size = resps)
-# plot(germany_geo) + plot(random_points, add=T)
-
-#### Draw respondents ####
-# Load population density data for Germany
-
-density_raster <- rast("Projects/ValuGaps/data/gis/BBSR_Landnutzung_Bevoelkerung/ewz250_2011_v1_multi.tif")
-
-
-##### Do plot of the pop density ####
-
-# density_df <- as.data.frame(density_raster, xy=TRUE)
-# 
-# density_df$cuts <- cut(density_df$ewz250_2011_v1_multi, breaks=c(0,5, 10, 15, 25,50,100,200,400, 800, 1200, 1600, 2000,4800),
-#                        include.lowest = TRUE) 
-# 
-# ggplot(data=density_df) +
-#   geom_raster(aes(x=x, y=y, fill=cuts)) +
-#   coord_equal() + 
-#   labs(fill="Population \ndensity") +
-#   xlab("") + ylab("") +
-#   scale_fill_viridis(discrete = TRUE, option="B", labels = c("5", "10", "15", "25", "50", "100", "200",
-#                                                              "400", "800", "1.200", "1.600", "3.200",
-#                                                              "4.800", "10.000", "> 10.000"))
-# 
-# ggsave("Figures/pop_density.png", width=5, height=4)
-  
-
-
-#### Draw cells from raster considering population density as probability weights ####
-
-random_points <- spatSample(density_raster, resps, method="weights", replace= TRUE, as.points=TRUE) 
-# set replace to TRUE to sample with replacement
-
-random_points_ls <- base::split(sample(random_points), 1:length(buffer_dist_ls))
-# Do a nice plot of the "respondents"
-
-# ggplot() +
-#   geom_sf(data=germany_geo) +
-#   geom_point(data=as.data.frame(random_points, geom="XY"), aes(x=x, y=y), size=0.5) +
-#   xlab("") + 
-#   ylab("")
-# 
-# ggsave("Figures/spat_dis_resp.png", width=5, height = 4, dpi="print")
-
-# Draw Buffer around random points 
-
-#### Split samples ####
-buffer_pts_ls <- list()
-
-buffer_pts_ls <- lapply(seq_along(random_points_ls), function(i) {
-  buffered <- buffer(random_points_ls[[i]], buffer_dist_ls[[i]])
-  buffered_sf <- st_as_sf(buffered, crs = crs(hnv))
-  buffered_sf
-})
-
-buffer_pts <- do.call(rbind, buffer_pts_ls)
-
-
-#buffer_pts <- buffer(random_points, buffer_dist)
-
-#buffer_pts <- st_as_sf(buffer_pts, crs=crs(hnv))
-
-### Equalize CRS's to avoid warning message; can be done here without problems as CRS is formally the same anyway
-
-crs(hnv) <- crs(corine_hnv) <- crs(protected_areas) <- crs(corine_pa)
-
-# Compute shares within Buffer with exact_extract and customized functions
-
-
-hnv_shares <- process_data(hnv, "HNV", buffer_pts_ls)
-
-hnv_possible <- process_data(corine_hnv, "HNV_possible", buffer_pts_ls)
-
-protected_shares <- process_data(protected_areas, "protected", buffer_pts_ls)
-
-protected_possible <- process_data(corine_pa, "Protected_possible", buffer_pts_ls)
-
-# Merge protected area shares and HNV shares
-land_cover_shares <- left_join(protected_shares, hnv_shares, by="ID")
-
-possible_shares <- left_join(protected_possible, hnv_possible, by="ID") 
-
-final_set <- left_join(land_cover_shares, possible_shares, by="ID") %>% 
-  mutate(HNV_on_corine = HNV/HNV_possible, Protected_on_corine = protected/Protected_possible,
-         buffer_dist_index = ceiling(row_number() / (resps / length(buffer_dist_ls))),
-         buffer_dist = map_dbl(buffer_dist_index, ~ buffer_dist_ls[[.]]),
-         SQ_HNV = (buffer_dist^2*HNV*pi)/10000,
-         SQ_Prot = (buffer_dist^2*protected*pi)/10000) 
-
-###################################################################################
-#### Overview of GIS shares ####
-
-# library(ggpubr)
-# pa_pl <- ggplot(data=final_set) +
-#   geom_density(aes(x=Protected_on_corine), fill="lightblue") +
-#   xlab("Protected areas (% of possible)") +
-#   xlim(0:1)
-# 
-# hnv_pl <- ggplot(data=final_set) +
-#   geom_density(aes(x=HNV_on_corine), fill="lightblue") +
-#   xlab("HNV (% of possible)") + xlim(0:1)
-# 
-# pa_pos <- ggplot(data=final_set) +
-#   geom_density(aes(x=Protected_possible), fill="yellow") +
-#   xlab("Protected possible share") + xlim(0:1)
-# 
-# hnv_pos <- ggplot(data=final_set) +
-#   geom_density(aes(x=HNV_possible), fill="yellow") +
-#   xlab("HNV possible share") + xlim(0:1)
-# 
-# pa_sq <- ggplot(data=final_set) +
-#   geom_density(aes(x=SQ_Prot), fill="lightgreen") +
-#   xlab("Protected areas SQ (ha)") + xlim(0,40000)
-# 
-# hnv_sq <- ggplot(data=final_set) +
-#   geom_density(aes(x=SQ_HNV), fill="lightgreen") +
-#   xlab("HNV SQ (ha)") + xlim(0,40000)
-# 
-# ggarrange(pa_pl, hnv_pl, pa_pos, hnv_pos, pa_sq, hnv_sq, nrow = 3, ncol=2)
-# 
-# ggsave("Figures/finalset_gis.png", dpi="print", width = 8, height = 6)
-###################################################################################
+# Functions to compute shares 
+
+sum_cover <- function(x){
+  list(x %>%
+         group_by(value) %>%
+         summarize(total_area = sum(coverage_area, na.rm=T)) %>%
+         mutate(proportion = total_area/sum(total_area[1:2])))
+  
+}
+
+process_data <- function(data, name_prefix, cover_data) {
+  result <- exact_extract(data, buffer_pts, fun=sum_cover, coverage_area=T, summarize_df=T)
+  col_name <- paste(name_prefix)
+  result <- as.data.frame(sapply(result, function(x) ifelse(length(x$proportion) >= 2, x$proportion[2], 0)))
+  colnames(result) <- col_name
+  result <- result %>% mutate(ID = 1:resps) %>% dplyr::select(ID, everything())
+  return(result)
+}
+
+
+##### Load in GIS data ####
+
+## Draw random points without incorporating actual pop density (old)
+# random_points <- st_sample(germany_geo_county, size = resps)
+# plot(germany_geo) + plot(random_points, add=T)
+
+#### Draw respondents ####
+# Load population density data for Germany
+
+density_raster <- rast("Projects/ValuGaps/data/gis/BBSR_Landnutzung_Bevoelkerung/ewz250_2011_v1_multi.tif")
+  
+
+
+#### Draw cells from raster considering population density as probability weights ####
+
+random_points <- spatSample(density_raster, resps, method="weights", replace= TRUE, as.points=TRUE) 
+# set replace to TRUE to sample with replacement
+
+random_points_ls <- base::split(sample(random_points), 1:length(buffer_dist_ls))
+
+# Draw Buffer around random points 
+
+#### Split samples ####
+buffer_pts_ls <- list()
+
+buffer_pts_ls <- lapply(seq_along(random_points_ls), function(i) {
+  buffered <- buffer(random_points_ls[[i]], buffer_dist_ls[[i]])
+  buffered_sf <- st_as_sf(buffered, crs = crs(hnv))
+  buffered_sf
+})
+
+buffer_pts <- do.call(rbind, buffer_pts_ls)
+
+
+
+### Equalize CRS's to avoid warning message; can be done here without problems as CRS is formally the same anyway
+
+crs(hnv) <- crs(corine_hnv) <- crs(protected_areas) <- crs(corine_pa)
+
+# Compute shares within Buffer with exact_extract and customized functions
+
+
+hnv_shares <- process_data(hnv, "HNV", buffer_pts_ls)
+
+hnv_possible <- process_data(corine_hnv, "HNV_possible", buffer_pts_ls)
+
+protected_shares <- process_data(protected_areas, "protected", buffer_pts_ls)
+
+protected_possible <- process_data(corine_pa, "Protected_possible", buffer_pts_ls)
+
+# Merge protected area shares and HNV shares
+land_cover_shares <- left_join(protected_shares, hnv_shares, by="ID")
+
+possible_shares <- left_join(protected_possible, hnv_possible, by="ID") 
+
+final_set <- left_join(land_cover_shares, possible_shares, by="ID") %>% 
+  mutate(HNV_on_corine = HNV/HNV_possible, Protected_on_corine = protected/Protected_possible,
+         buffer_dist_index = ceiling(row_number() / (resps / length(buffer_dist_ls))),
+         buffer_dist = map_dbl(buffer_dist_index, ~ buffer_dist_ls[[.]]),
+         SQ_HNV = (buffer_dist^2*HNV*pi)/10000,
+         SQ_Prot = (buffer_dist^2*protected*pi)/10000) 
+
+###################################################################################
+#### Overview of GIS shares ####
+
+# library(ggpubr)
+# pa_pl <- ggplot(data=final_set) +
+#   geom_density(aes(x=Protected_on_corine), fill="lightblue") +
+#   xlab("Protected areas (% of possible)") +
+#   xlim(0:1)
+# 
+# hnv_pl <- ggplot(data=final_set) +
+#   geom_density(aes(x=HNV_on_corine), fill="lightblue") +
+#   xlab("HNV (% of possible)") + xlim(0:1)
+# 
+# pa_pos <- ggplot(data=final_set) +
+#   geom_density(aes(x=Protected_possible), fill="yellow") +
+#   xlab("Protected possible share") + xlim(0:1)
+# 
+# hnv_pos <- ggplot(data=final_set) +
+#   geom_density(aes(x=HNV_possible), fill="yellow") +
+#   xlab("HNV possible share") + xlim(0:1)
+# 
+# pa_sq <- ggplot(data=final_set) +
+#   geom_density(aes(x=SQ_Prot), fill="lightgreen") +
+#   xlab("Protected areas SQ (ha)") + xlim(0,40000)
+# 
+# hnv_sq <- ggplot(data=final_set) +
+#   geom_density(aes(x=SQ_HNV), fill="lightgreen") +
+#   xlab("HNV SQ (ha)") + xlim(0,40000)
+# 
+# ggarrange(pa_pl, hnv_pl, pa_pos, hnv_pos, pa_sq, hnv_sq, nrow = 3, ncol=2)
+# 
+# ggsave("Figures/finalset_gis.png", dpi="print", width = 8, height = 6)
+###################################################################################

Projects/ValuGaps/designsyntax/design syntax.ngs

@@ -140,6 +140,7 @@ $
 ?Efficient Design
 
 
+
 design
 ;alts = alt1*, alt2*, alt3
 ;rows = 36
@@ -149,60 +150,11 @@ design
 
 ;con
 ;model:
-U(alt1) = b1[-1.2] + b2[0.1] * B[0,1] + b3[0.4] * C[0,1]  + b4[0.3] * D[0,1] + b5[0.02] * P[35:83:6]/
-U(alt2) =  b6[-1.4] + b2               * B      + b3            * C       + b4            * D      + b5                       * P
+U(alt1) = b1[-1.2] + bPROT[0.1] * B[0,10,20] + bHNV[0.1] * C[0,10,20]  + bPROTACC[0.2] * D[0,1,2] + bHNVACC[0.1] * E[0,1] +bPRICE[-0.02] *P[40,80,120,160,200,300]/ 
+U(alt2) = b1       + bPROT      * B          + bHNV      * C           + bPROTACC      * D        + bHNVACC      * E      +bPRICE        *P[40,80,120,160,200,300] 
 ?U(alt3) = 
 
-;formatTitle = 'Scenario <scenarionumber>'
-;formatTableDimensions = 3, 6
-;formatTable:
-1,1 = '' /
-1,2 = 'Results' /
-1,3 = 'Beratung' /
-1,4 = 'Partner' /
-1,5 = 'Kompensation' /
-1,6 = 'Choice question&:' /
-2,1 = 'alt1' /
-2,2 = '<alt1.b>' /
-2,3 = '<alt1.c>' /
-2,4 = '<alt1.d>' /
-2,5 = '<alt1.p>' /
-2,6 = '' /
-3,1 = 'alt2' /
-3,2 = '<alt2.b>' /
-3,3 = '<alt2.c>' /
-3,4 = '<alt2.d>' /
-3,5 = '<alt2.p>' /
-3,6 = ''
-;formatTableStyle:
-1,1 = 'default' /
-1,2 = 'headingattribute' /
-1,3 = 'headingattribute' /
-1,4 = 'headingattribute' /
-1,5 = 'headingattribute' /
-1,6 = 'headingattribute' /
-2,1 = 'heading1' /
-2,2 = 'body1' /
-2,3 = 'body1' /
-2,4 = 'body1' /
-2,5 = 'body1' /
-2,6 = 'choice1' /
-3,1 = 'heading2' /
-3,2 = 'body2' /
-3,3 = 'body2' /
-3,4 = 'body2' /
-3,5 = 'body2' /
-3,6 = 'choice2'
-;formatStyleSheet = Default.css
-;formatAttributes:
-alt1.b(0=Results, 1=Action) /
-alt1.c(0=Keine Beratung, 1=Mit Beratung) /
-alt1.d(0=Keine Partner, 1=Mit Partner) /
-alt1.p(35=#0, 41=#0, 47=#0, 53=#0, 59=#0, 65=#0, 71=#0, 77=#0, 83=#0) /
-alt2.b(0=Results, 1=Action) /
-alt2.c(0=Keine Beratung, 1=Mit Beratung) /
-alt2.d(0=Keine Partner, 1=Mit Partner) /
-alt2.p(35=#0, 41=#0, 47=#0, 53=#0, 59=#0, 65=#0, 71=#0, 77=#0, 83=#0)
+
 $
 
 

Projects/ValuGaps/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))
-####
-
-##########################################################
-# Check how many cells are both, HNV and protected areas #
-
-# pa_resample <- resample(protected_areas, hnv) # bring both rasters to same extent
-# pa_resample <- ifel(pa_resample==1, 3, 0) # recode such that 2 is equal to overlap
-# test_overlap <- pa_resample - hnv
-# test_overlap <- ifel(test_overlap==2, 1, 0)  #recode overlap cells at 1, rest 0
-# plot(test_overlap)
-# summary(test_overlap)
-# mean(values(test_overlap), na.rm=T) # 0.7% of Germany's landcover are overlapping cells
-# 
-# ggplot(data=as.data.frame(test_overlap==1, xy=TRUE)) +
-#   geom_raster(aes(x=x, y=y, fill=value)) +
-#   scale_fill_manual(values=c("grey80", "darkred"), name="Overlap") +
-#   coord_equal()
-# ggsave("Figures/overlap.png", width=10, height = 8)
-##########################################################
-
-
-#### Define simulation parameters ####
-
-resps = 96  # number of respondents 
-nosim=3  # 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(
-  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 )
-
-# u<-list(
-# v1 =V.1~ basc + bb*alt1.b+ bc * alt1.c + bp * alt1.p , 
-# v2 =V.2~ basc + bb*alt2.b+ bc * alt2.c + bp * alt2.p ,
-# v3 =V.3~ 0)
-
-
-## this function will be called before the data is simulated
-sou_gis <- function() {
-  
-  source("Projects/ValuGaps/code/GIS_data.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=3  # 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(
+  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 )
+
+# u<-list(
+# v1 =V.1~ basc + bb*alt1.b+ bc * alt1.c + bp * alt1.p , 
+# v2 =V.2~ basc + bb*alt2.b+ bc * alt2.c + bp * alt2.p ,
+# v3 =V.3~ 0)
+
+
+## this function will be called before the data is simulated
+sou_gis <- function() {
+  
+  source("Projects/ValuGaps/code/GIS_data.R")
+  
+}

functions.R

@@ -15,6 +15,7 @@ sim_choice <- function(designfile, no_sim=10, respondents=330, mnl_U,utils=u ) {
   require("purrr")
   require("ggplot2")
   require("formula.tools")  
+  require("rlang")
    
   
   mnl_U <-paste(map_chr(u,as.character,keep.source.attr = TRUE),collapse = "",";") %>%

simulationcore_purrr.R

@@ -6,9 +6,11 @@ library("formula.tools")
 library(ggplot2)
 library(stringr)
 
-source("Projects/SE_AGRI/parameters_SE Design-Agri.R")
+
+#source("Projects/SE_DRIVE/parameters_SE_DRIVE.R")
+#source("Projects/SE_AGRI/parameters_SE Design-Agri.R")
 #source("Projects/ValuGaps/parameters_valugaps.R",echo = TRUE)
-#source(params$file)
+source(params$file)
 
 
 designfile<-list.files(designpath,full.names = T)