From 1d0fca1fd3e36e6cbc50efcaa728fa61ff3d4a14 Mon Sep 17 00:00:00 2001
From: Julian Sagebiel <julian.sagebiel@idiv.de>
Date: Fri, 22 Sep 2023 17:21:54 +0200
Subject: [PATCH] working version with small changes
---
Projects/ValuGaps/code/GIS_data.R | 261 ++++++++--------
.../ValuGaps/designsyntax/design syntax.ngs | 56 +---
Projects/ValuGaps/parameters_valugaps.R | 278 +++++++++---------
functions.R | 1 +
simulationcore_purrr.R | 6 +-
5 files changed, 256 insertions(+), 346 deletions(-)
diff --git a/Projects/ValuGaps/code/GIS_data.R b/Projects/ValuGaps/code/GIS_data.R
index ac73130..12e44f9 100644
--- a/Projects/ValuGaps/code/GIS_data.R
+++ b/Projects/ValuGaps/code/GIS_data.R
@@ -1,146 +1,115 @@
-# 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)
+###################################################################################
diff --git a/Projects/ValuGaps/designsyntax/design syntax.ngs b/Projects/ValuGaps/designsyntax/design syntax.ngs
index a15be51..e64339c 100644
--- a/Projects/ValuGaps/designsyntax/design syntax.ngs
+++ b/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)
+
$
diff --git a/Projects/ValuGaps/parameters_valugaps.R b/Projects/ValuGaps/parameters_valugaps.R
index 9df3ee3..0a48b33 100644
--- a/Projects/ValuGaps/parameters_valugaps.R
+++ b/Projects/ValuGaps/parameters_valugaps.R
@@ -1,146 +1,132 @@
-
-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")
+
+}
diff --git a/functions.R b/functions.R
index 1feaa7e..1a3bbe9 100644
--- a/functions.R
+++ b/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 = "",";") %>%
diff --git a/simulationcore_purrr.R b/simulationcore_purrr.R
index 67f02e0..345b601 100644
--- a/simulationcore_purrr.R
+++ b/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)
--
GitLab