From bfffdae7194b0986ceea3b332de903388cfbdf3e Mon Sep 17 00:00:00 2001
From: dj44vuri <julian.sagebiel@idiv.de>
Date: Mon, 11 Dec 2023 14:25:12 +0100
Subject: [PATCH] add valugaps without data
---
.gitignore | 2 +-
Projects/ValuGaps/Designs/ort2attr.ngd | 156 ++++++
Projects/ValuGaps/code/BFN_data.R | 53 ++
Projects/ValuGaps/code/GIS_data.R | 115 +++++
Projects/ValuGaps/code/makes_things_faster.R | 24 +
Projects/ValuGaps/code/maps.R | 469 ++++++++++++++++++
.../ValuGaps/designsyntax/design syntax.ngs | 234 +++++++++
Projects/ValuGaps/parameters_valugaps.R | 129 +++++
8 files changed, 1181 insertions(+), 1 deletion(-)
create mode 100644 Projects/ValuGaps/Designs/ort2attr.ngd
create mode 100644 Projects/ValuGaps/code/BFN_data.R
create mode 100644 Projects/ValuGaps/code/GIS_data.R
create mode 100644 Projects/ValuGaps/code/makes_things_faster.R
create mode 100644 Projects/ValuGaps/code/maps.R
create mode 100644 Projects/ValuGaps/designsyntax/design syntax.ngs
create mode 100644 Projects/ValuGaps/parameters_valugaps.R
diff --git a/.gitignore b/.gitignore
index 772736c..b3b73fd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,5 +6,5 @@ dce_simulation_tool.Rproj
modeloutput/
*.html
simulation_output_files/
-Projects/ValuGaps/
+Projects/ValuGaps/data
Projects/feedadditives/feedadd_alldestest_files
\ No newline at end of file
diff --git a/Projects/ValuGaps/Designs/ort2attr.ngd b/Projects/ValuGaps/Designs/ort2attr.ngd
new file mode 100644
index 0000000..56a428b
--- /dev/null
+++ b/Projects/ValuGaps/Designs/ort2attr.ngd
@@ -0,0 +1,156 @@
+Design Choice situation alt1.b alt1.c alt1.p alt2.b alt2.c alt2.p Block
+1 1 2 3 40 3 2 120 5
+1 2 3 2 80 2 3 120 5
+1 3 0 1 120 1 0 120 5
+1 4 1 0 160 0 1 120 5
+1 5 2 0 120 3 1 80 2
+1 6 3 1 160 2 0 80 2
+1 7 0 2 40 1 3 80 2
+1 8 1 3 80 0 2 80 2
+1 9 1 2 120 0 3 40 6
+1 10 0 3 160 1 2 40 6
+1 11 3 0 40 2 1 40 6
+1 12 2 1 80 3 0 40 6
+1 13 3 3 160 3 3 160 1
+1 14 2 2 120 2 2 160 1
+1 15 1 1 80 1 1 160 1
+1 16 0 0 40 0 0 160 1
+1 17 3 3 80 1 1 80 3
+1 18 2 2 40 0 0 80 3
+1 19 1 1 160 3 3 80 3
+1 20 0 0 120 2 2 80 3
+1 21 2 2 80 1 1 160 5
+1 22 3 3 40 0 0 160 5
+1 23 0 0 160 3 3 160 5
+1 24 1 1 120 2 2 160 5
+1 25 0 2 80 3 0 40 1
+1 26 1 3 40 2 1 40 1
+1 27 2 0 160 1 2 40 1
+1 28 3 1 120 0 3 40 1
+1 29 0 1 40 2 1 160 6
+1 30 1 0 80 3 0 160 6
+1 31 2 3 120 0 3 160 6
+1 32 3 2 160 1 2 160 6
+1 33 0 2 160 0 1 120 4
+1 34 1 3 120 1 0 120 4
+1 35 2 0 80 2 3 120 4
+1 36 3 1 40 3 2 120 4
+1 37 0 3 80 2 3 80 4
+1 38 1 2 40 3 2 80 4
+1 39 2 1 160 0 1 80 4
+1 40 3 0 120 1 0 80 4
+1 41 2 3 160 2 0 40 2
+1 42 3 2 120 3 1 40 2
+1 43 0 1 80 0 2 40 2
+1 44 1 0 40 1 3 40 2
+1 45 1 2 160 2 0 120 3
+1 46 0 3 120 3 1 120 3
+1 47 3 0 80 0 2 120 3
+1 48 2 1 40 1 3 120 3
+1 49 2 1 40 1 3 40 3
+1 50 3 0 80 0 2 40 3
+1 51 0 3 120 3 1 40 3
+1 52 1 2 160 2 0 40 3
+1 53 1 0 40 1 3 120 2
+1 54 0 1 80 0 2 120 2
+1 55 3 2 120 3 1 120 2
+1 56 2 3 160 2 0 120 2
+1 57 3 0 120 1 0 160 4
+1 58 2 1 160 0 1 160 4
+1 59 1 2 40 3 2 160 4
+1 60 0 3 80 2 3 160 4
+1 61 3 1 40 3 2 40 4
+1 62 2 0 80 2 3 40 4
+1 63 1 3 120 1 0 40 4
+1 64 0 2 160 0 1 40 4
+1 65 3 2 160 1 2 80 6
+1 66 2 3 120 0 3 80 6
+1 67 1 0 80 3 0 80 6
+1 68 0 1 40 2 1 80 6
+1 69 3 1 120 0 3 120 1
+1 70 2 0 160 1 2 120 1
+1 71 1 3 40 2 1 120 1
+1 72 0 2 80 3 0 120 1
+1 73 1 1 120 2 2 80 5
+1 74 0 0 160 3 3 80 5
+1 75 3 3 40 0 0 80 5
+1 76 2 2 80 1 1 80 5
+1 77 0 0 120 2 2 160 3
+1 78 1 1 160 3 3 160 3
+1 79 2 2 40 0 0 160 3
+1 80 3 3 80 1 1 160 3
+1 81 0 0 40 0 0 80 1
+1 82 1 1 80 1 1 80 1
+1 83 2 2 120 2 2 80 1
+1 84 3 3 160 3 3 80 1
+1 85 2 1 80 3 0 120 6
+1 86 3 0 40 2 1 120 6
+1 87 0 3 160 1 2 120 6
+1 88 1 2 120 0 3 120 6
+1 89 1 3 80 0 2 160 2
+1 90 0 2 40 1 3 160 2
+1 91 3 1 160 2 0 160 2
+1 92 2 0 120 3 1 160 2
+1 93 1 0 160 0 1 40 5
+1 94 0 1 120 1 0 40 5
+1 95 3 2 80 2 3 40 5
+1 96 2 3 40 3 2 40 5
+||||||||||
+design
+;alts = alt1*, alt2*, alt3
+;rows = 72
+;block = 6
+;orth = sim
+
+
+;con
+;model:
+U(alt1) = b1 + b2 * B[0,1,2,3] + b3 * C[0,1,2,3] + b5 * P[40,80,120,160]/
+U(alt2) = b1 + b2 * B + b3 * C + b5 * P
+
+
+
+
+;formatTitle = 'Scenario <scenarionumber> Block <blocknumber>'
+;formatTableDimensions = 3, 5
+;formatTable:
+1,1 = '' /
+1,2 = 'Share of protected areas' /
+1,3 = 'Share of HNV' /
+1,4 = 'Cost per year' /
+1,5 = 'Choice question&:' /
+2,1 = 'alt1' /
+2,2 = '<alt1.b>' /
+2,3 = '<alt1.c>' /
+2,4 = '<alt1.p>' /
+2,5 = '' /
+3,1 = 'alt2' /
+3,2 = '<alt2.b>' /
+3,3 = '<alt2.c>' /
+3,4 = '<alt2.p>' /
+3,5 = ''
+;formatTableStyle:
+1,1 = 'default' /
+1,2 = 'headingattribute' /
+1,3 = 'headingattribute' /
+1,4 = 'headingattribute' /
+1,5 = 'headingattribute' /
+2,1 = 'heading1' /
+2,2 = 'body1' /
+2,3 = 'body1' /
+2,4 = 'body1' /
+2,5 = 'choice1' /
+3,1 = 'heading2' /
+3,2 = 'body2' /
+3,3 = 'body2' /
+3,4 = 'body2' /
+3,5 = 'choice2'
+;formatStyleSheet = Default.css
+;formatAttributes:
+alt1.b(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt1.c(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt1.p(40=#, 80=#, 120=#, 160=#) /
+alt2.b(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt2.c(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt2.p(40=#, 80=#, 120=#, 160=#)
+$
\ No newline at end of file
diff --git a/Projects/ValuGaps/code/BFN_data.R b/Projects/ValuGaps/code/BFN_data.R
new file mode 100644
index 0000000..599f0c3
--- /dev/null
+++ b/Projects/ValuGaps/code/BFN_data.R
@@ -0,0 +1,53 @@
+# Only needed to be run one time to get bfn data
+# can be adjusted to load data for other layers (nationalparks etc.)
+library(xml2)
+library(terra)
+library(sf)
+library(tidyverse)
+library(gdalUtilities)
+
+url <- "https://geodienste.bfn.de/ogc/wfs/schutzgebiet?REQUEST=GetCapabilities&SERVICE=WFS&VERSION=2.0.0"
+wms_capabilities <- read_xml(url)
+
+layers <- st_layers(url)
+
+data_nasc <- read_sf(url, layer= "bfn_sch_Schutzgebiet:Naturschutzgebiete")
+
+data_nasc <- data_nasc %>% mutate(value=1)
+
+ggplot(data=data_nasc) +
+ #geom_sf(data=germany_geo) +
+ geom_sf(fill="darkgreen", col="darkgreen")
+
+#ggsave("Figures/bfn_map.png", width = 5, height = 4)
+
+
+# Rasterize data #
+
+# Issue with conversion use function to convert data first before vectorizing
+ensure_multipolygons <- function(X) {
+ tmp1 <- tempfile(fileext = ".gpkg")
+ tmp2 <- tempfile(fileext = ".gpkg")
+ st_write(X, tmp1)
+ ogr2ogr(tmp1, tmp2, f = "GPKG", nlt = "MULTIPOLYGON")
+ Y <- st_read(tmp2)
+ st_sf(st_drop_geometry(X), geom = st_geometry(Y))
+}
+
+data_nasc_con <- ensure_multipolygons(data_nasc)
+
+data_nasc_con <- st_transform(data_nasc_con, crs="epsg:3035")
+
+vec_nasc <- vect(data_nasc_con)
+
+empty_raster <- rast(vec_nasc, resolution(1,1))
+res(empty_raster) <- c(100, 100)
+
+raster_nasc <- rasterize(vec_nasc, empty_raster, "value")
+
+plot(raster_nasc, col="darkgreen")
+
+saveRDS(raster_nasc, "data/gis/protected_bfn_raster.rds")
+
+# Check if conversion and storage worked #
+raster_bfn <- readRDS("data/gis/protected_bfn_raster.rds") # works
diff --git a/Projects/ValuGaps/code/GIS_data.R b/Projects/ValuGaps/code/GIS_data.R
new file mode 100644
index 0000000..12e44f9
--- /dev/null
+++ b/Projects/ValuGaps/code/GIS_data.R
@@ -0,0 +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")
+
+
+
+#### 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/code/makes_things_faster.R b/Projects/ValuGaps/code/makes_things_faster.R
new file mode 100644
index 0000000..e591a2a
--- /dev/null
+++ b/Projects/ValuGaps/code/makes_things_faster.R
@@ -0,0 +1,24 @@
+# Do some time testing
+tic()
+crs(buffer_pts) <- crs(hnv)
+hnv_shares <- terra::extract(hnv, buffer_pts, fun= mean, na.rm =T, weights = T)
+toc()
+
+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])))
+
+}
+
+# More fast but not correct for points at the edge
+tic()
+poly <- st_as_sf(buffer_pts, crs=crs(hnv))
+test_share <- exact_extract(hnv, poly, fun=sum_cover, coverage_area=T, summarize_df=T)
+test_share <- as.data.frame(sapply(test_share, function(x) ifelse(length(x$proportion) >= 2, x$proportion[2], 0)))
+colnames(test_share) <- "HNV"
+test_share <- test_share %>% mutate(ID = 1:resps)
+toc()
+
+check <- (hnv_shares$eea_r_3035_100_m_hnv.farmland.ac_2012_v1_r00 - test_share$`sapply(test_share, function(x) x$proportion[2])`)
diff --git a/Projects/ValuGaps/code/maps.R b/Projects/ValuGaps/code/maps.R
new file mode 100644
index 0000000..f2b7b70
--- /dev/null
+++ b/Projects/ValuGaps/code/maps.R
@@ -0,0 +1,469 @@
+library(tidyverse)
+library(tidylog)
+library(apollo)
+library(sf)
+library(terra)
+library(exactextractr)
+library(ggspatial)
+library(viridis)
+
+#### 1 Define parameters and functions ####
+
+buffer_dist = 15000 # distance we use in dce
+
+### Summarizing function
+
+sum_cover <- function(x){
+ list(x %>%
+ group_by(value) %>%
+ summarize(total_area = sum(coverage_area)) %>%
+ mutate(proportion = total_area/sum(total_area)))
+
+}
+
+#### 2 Load in model data ####
+map_data <- readRDS("output/96_3runs_4designs_mixl_2023-08-14.RDS")
+
+# Use estimates from first simulated model for example map
+
+model1 <- as.data.frame(map_data$ort2attr[[1]]$estimate)
+
+model1 <- rownames_to_column(model1)
+colnames(model1) <- c("Parameter", "Estimate")
+
+#### 3 Read in data for Germany with administrative borders ####
+
+# https://gadm.org/download_country.html #
+germany_county <- read_sf("data/gadm41_DEU_shp", "gadm41_DEU_2") # county borders
+germany_gemeinde <- read_sf("data/gadm41_DEU_shp", "gadm41_DEU_3") # Gemeinde borders
+
+
+#### 4 Read in CORINE Land cover data ####
+
+corine <- rast("data/gis/corine_raster_2018/DATA/U2018_CLC2018_V2020_20u1.tif")
+germany_geo <- st_transform(germany_county$geometry, crs=crs(corine)) # used to crop corine data to Germany
+corine <- crop(corine, germany_geo)
+corine <- mask(corine, vect(germany_geo))
+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))
+
+#### 5 High Nature Value (HNV) ####
+
+# Calculate share of HNV for each county
+hnv <- rast("data/gis/hnv_germany.tif")
+
+###### Overlap of HNV and hnv_corine ####
+
+hnv_true <- hnv > 0
+hnv_filtered <- mask(hnv, hnv_true, maskvalue=0)
+
+hnv_on_corine <- intersect(corine_hnv, hnv_filtered)
+hnv_on_corine <- mask(hnv_on_corine, vect(germany_geo)) #hnv only on the categories where we would expect it
+
+corine_hnv_true <- corine_hnv > 0
+corine_hnv_test <- mask(corine_hnv, corine_hnv_true, maskvalue=0)
+
+germany_county <- st_transform(germany_county, crs(corine_hnv))
+germany_gemeinde <- st_transform(germany_gemeinde, crs(corine_hnv))
+
+##### Calculate actual and possible HNV share for all counties ######
+
+counties <- germany_county$NAME_2
+hnv_list <- list()
+hnv_possible_list <- list()
+
+z= 1
+
+
+for (county in counties) {
+ temp_county <- germany_county %>% filter(NAME_2 == counties[z])
+ hnv_possible_temp <- as.data.frame(
+ exact_extract(corine_hnv, temp_county, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)) # caclulate possible area
+ hnv_possible_list[[county]] <- hnv_possible_temp
+ hnv_temp <- as.data.frame(
+ exact_extract(hnv, temp_county, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)) # calculate HNV share on total area
+ hnv_list[[county]] <- hnv_temp
+ print(paste0(n_distinct(germany_county$NAME_2) - z, " remaining"))
+ z = z+1
+}
+
+hnv_shares_counties <- as.data.frame(sapply(hnv_list, function(x) x$proportion[2]))
+hnv_shares_counties <- rownames_to_column(hnv_shares_counties)
+hnv_shares_counties$possible <- as.numeric(sapply(hnv_possible_list, function(x) x$proportion[2]))
+colnames(hnv_shares_counties) <- c("County", "HNV_share", "HNV_Possible")
+hnv_shares_counties <- hnv_shares_counties %>% mutate(HNV_share = replace_na(HNV_share, 0),
+ HNV_on_corine = HNV_share/HNV_Possible)
+
+
+##### check total HNV area in relation to Germany and Germany agricultural area ####
+hnv_area <- as.data.frame(sapply(hnv_list, function(x) x$total_area[2]))
+names(hnv_area) <- "Area"
+hnv_area_sum <- hnv_area %>% summarize(Total_area = sum(Area, na.rm=T)/10000)
+hnv_area_possible <- as.data.frame(sapply(hnv_possible_list, function(x) x$total_area[2]))
+names(hnv_area_possible) <- "Area"
+hnv_area_sum_possible <- hnv_area_possible %>% summarize(Total_area = sum(Area, na.rm=T)/10000)
+
+share_germany_hnv <- as.numeric(hnv_area_sum)/as.numeric(hnv_area_sum_possible)
+
+
+
+###### Calculate WTP based on share for HNV ####
+
+wtp_hnv <- function(x) {-(model1$Estimate[[4]] + 2*model1$Estimate[5]*x)/model1$Estimate[6]}
+
+hnv_shares_counties <- hnv_shares_counties %>% mutate(WTP_HNV = wtp_hnv(HNV_share*100))
+
+###### Create HNV map ####
+
+# Merge county shapes with data frame
+
+germany_county <- germany_county %>% rename("County" = NAME_2)
+merge_hnv <- hnv_shares_counties %>% left_join(germany_county, by="County") %>%
+ select(County, HNV_share, WTP_HNV, geometry)
+
+
+ggplot(data=germany_county) +
+ geom_sf(fill="gray98", alpha=1) +
+ geom_sf(data=merge_hnv, aes(fill=WTP_HNV, geometry=geometry), size=2) +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ labs(fill= "WTP (€)") +
+ scale_fill_gradientn(colours = c("cornsilk1", "royalblue3")) +
+ theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm")) +
+ ggtitle("WTP for a 1% increase in HNV")
+
+ggsave("Figures/hnv_cont.png", width=6, height=5, dpi="print")
+
+## discontinous map ##
+
+# quantile(merge_hnv$WTP_HNV)
+# merge_hnv$cuts <- cut(merge_hnv$WTP_HNV, breaks = c(14, 22, 30, 38, 47))
+#
+#
+# ggplot(data=germany_county) +
+# geom_sf(fill="gray98") +
+# geom_sf(data=merge_hnv, aes(fill=cuts, geometry=geometry), size=2) +
+# annotation_scale(location = "bl", width_hint = 0.12,plot_unit="m", line_width = 0.5) +
+# annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+# style = north_arrow_fancy_orienteering) +
+# labs(fill= "WTP (€)") +
+# scale_fill_manual(values = c("red3","orange", "gold1", "chartreuse3", "forestgreen"),
+# labels = c("< -17", " -17 - 0", "0 - 23", "23 - 46", "46 - 64"))
+#
+# ggsave("Figures/hnv_disc.png", dpi="print", width=7, height=5)
+
+
+#### 6 Protected Areas ####
+
+#### Protected Areas and CORINE
+
+protected_area_categories <- c(10, 12:41) # remove artifical 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))
+
+##### Load PA data ######
+
+protected_areas_all <- readRDS("data/gis/germany_protected.rds")
+
+categories_to_drop <- c(0,1,2,3, 4, 6, 8, 11) # define categories that should not be included in protected areas
+
+protected_areas <- protected_areas_all
+protected_areas[protected_areas %in% categories_to_drop] <-NA # Make binary categories
+protected_areas <- ifel(is.na(protected_areas), 0, 1)
+protected_areas <- mask(protected_areas, vect(germany_geo))
+
+prot_list <- list()
+prot_possible_list <- list()
+
+z= 1
+for (county in counties) {
+ temp_county <- germany_county %>% filter(County == counties[z])
+ prot_corine_temp <- as.data.frame(
+ exact_extract(corine_pa, temp_county, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover))
+ prot_possible_list[[county]] <- prot_corine_temp
+ prot_temp <- as.data.frame(
+ exact_extract(protected_areas, temp_county, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover))
+ prot_list[[county]] <- prot_temp
+ z = z+1
+}
+prot_counties <- as.data.frame(sapply(prot_list, function(x) x$proportion[2]))
+prot_counties <- rownames_to_column(prot_counties)
+prot_counties$Prot_Possible <- as.numeric(sapply(prot_possible_list, function(x) x$proportion[2]))
+colnames(prot_counties) <- c("County", "Prot_share", "Prot_possible")
+prot_counties <- prot_counties %>% mutate(Prot_share = replace_na(Prot_share , 0))
+prot_counties <- prot_counties %>% mutate(Prot_share_on_corine = Prot_share/Prot_possible)
+
+
+###### Calculate WTP based on share for Protected areas ####
+
+wtp_prot <- function(x) {-(model1$Estimate[[2]] + 2*model1$Estimate[3]*x)/model1$Estimate[6]}
+
+prot_counties <- prot_counties %>% mutate(WTP_Prot = wtp_prot(Prot_share*100))
+
+merge_prot <- prot_counties %>% left_join(germany_county, by="County") %>%
+ select(County, Prot_share, WTP_Prot, geometry)
+
+ggplot(data=germany_county) +
+ geom_sf(fill="gray98", alpha=1) +
+ geom_sf(data=merge_prot, aes(fill=WTP_Prot, geometry=geometry), size=2) +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ labs(fill= "WTP (€)") +
+ #scale_fill_viridis(option="C")
+ scale_fill_gradientn(colours = c("cornsilk1", "darkgreen")) +
+ theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"))
+
+
+ggsave("Figures/prot_cont.png",width=6, height=5, dpi="print")
+
+# make map for Gemeinde protected
+
+gemeinden <- germany_gemeinde$NAME_3
+gemeinden_code <- germany_gemeinde$CC_3 # use unique code since names are doubling for some gemeinden
+prot_list_g <- list()
+
+z= 1
+for (gemeinde in gemeinden) {
+ temp_county <- germany_gemeinde %>% filter(CC_3 == gemeinden_code[z])
+ prot_temp <- as.data.frame(
+ exact_extract(protected_areas, temp_county, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover))
+ prot_list_g[[gemeinde]] <- prot_temp
+ z = z+1
+}
+
+prot_gemeinde <- as.data.frame(sapply(prot_list_g, function(x) x$proportion[2]))
+prot_gemeinde <- rownames_to_column(prot_gemeinde)
+colnames(prot_gemeinde) <- c("Gemeinde", "Prot_share")
+prot_gemeinde <- prot_gemeinde %>% mutate(Prot_share = replace_na(Prot_share , 0))
+
+
+prot_gemeinde <- prot_gemeinde %>% mutate(WTP_Prot = wtp_prot(Prot_share*100))
+
+germany_gemeinde <- germany_gemeinde %>% rename("Gemeinde" = NAME_3)
+
+merge_prot_g <- prot_gemeinde %>% left_join(germany_gemeinde, by="Gemeinde") %>%
+ select(Gemeinde, Prot_share, WTP_Prot, geometry)
+
+ggplot(data=germany_gemeinde) +
+ geom_sf(fill="gray98", alpha=1) +
+ geom_sf(data=merge_prot_g, aes(fill=WTP_Prot, geometry=geometry), size=2, col="NA") +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ labs(fill= "WTP (€)") +
+ #scale_fill_viridis(option="C")
+ scale_fill_gradientn(colours = c("indianred", "white", "darkseagreen", "darkgreen")) +
+ ggtitle("WTP for a 1% increase in PA Share") +
+ theme(plot.margin = unit(c(0.1,0.1,0.1,0.1), "cm"))
+
+ggsave("Figures/prot_g_cont.png", dpi="print", width=6, height=5)
+
+
+### Do plot for protected shares gemeinde ###
+
+ggplot(data=merge_prot_g) +
+ geom_density(aes(x=Prot_share*100), col="darkseagreen", fill="darkseagreen", alpha=0.5) +
+ xlab("Share of Protected Areas (%)") +
+ ggtitle("Overview Protected Area Shares Gemeinden")
+
+ggsave("Figures/prot_share_g.png", width = 7, height = 5, dpi="print")
+
+
+# Do plot for HNV Share counties
+ggplot(data=merge_hnv) +
+ geom_density(aes(x=HNV_share*100), col="royalblue", fill="royalblue", alpha=0.5) +
+ xlab("HNV Share (%)") +
+ ggtitle("Overview HNV Shares Counties")
+
+ggsave("Figures/hnv_share.png", width = 7, height = 5, dpi="print")
+
+
+### 7 Raster values WTP #####
+
+##### Scaling ####
+
+scale_factor <- pi*(buffer_dist/1000)^2 # Calculate total area that is valued in the DCE as scaling factor
+
+###### Raster HNV (need to work on this further) ####
+#inefficient as hell
+
+agg_hnv <- aggregate(hnv, fact=10, fun=mean)
+
+check_list <- c(unlist(values(agg_hnv)))
+
+
+# try something, check this
+agg_hnv_wtp <- -(model1$Estimate[[4]] + 2*model1$Estimate[5]*(check_list*100))/model1$Estimate[6]
+
+values(agg_hnv) <- agg_hnv_wtp/scale_factor
+
+
+
+agg_hnv_gg <- as.data.frame(agg_hnv, xy=TRUE)
+colnames(agg_hnv_gg) <- c("x", "y", "WTP")
+
+# somehow raster only works sometimes, tile works but looks not as nice
+ggplot(data=agg_hnv_gg) +
+ geom_tile(aes(x = x, y = y, fill = WTP)) +
+ coord_quickmap() +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ #scale_fill_viridis(option="C", direction = 1)
+ scale_fill_gradientn(colours = c("gold", "cornsilk1", "royalblue3"))
+
+ggsave("Figures/hnv_wtp_raster.png", width=6, height=5, dpi="print")
+
+###### Raster protected ####
+
+agg_prot <- aggregate(protected_areas, fact=10, fun=mean)
+
+agg_prot <- mask(agg_prot, germany_county)
+
+values_vec <- c(unlist(values(agg_prot)))
+
+prot_wtp <- -(model1$Estimate[[2]] + 2*model1$Estimate[3]*values_vec*100)/model1$Estimate[6]
+
+values(agg_prot) <- prot_wtp
+
+
+agg_prot_gg <- as.data.frame(agg_prot, xy=TRUE)
+colnames(agg_prot_gg) <- c("x", "y", "WTP")
+
+# Plot protected
+ggplot(data=agg_prot_gg) +
+ geom_raster(aes(x = x, y = y, fill = WTP)) +
+ coord_equal() +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ #scale_fill_viridis(option="C", direction = 1)
+ scale_fill_gradientn(colours = c("cornsilk", "white", "darkseagreen", "darkgreen"))
+
+ggsave("Figures/prot_wtp_raster.png", width=6, height=5, dpi="print")
+
+#### 8 Load in Pop density data ####
+
+density_raster <- rast("data/gis/BBSR_Landnutzung_Bevoelkerung/ewz250_2011_v1_multi.tif")
+
+dens_rast_1000 <- aggregate(density_raster, fact=4, fun=sum)
+
+dens_gg <- as.data.frame(dens_rast_1000, xy=TRUE)
+colnames(dens_gg) <- c("x", "y", "Density")
+
+
+
+ggplot(data=dens_gg) +
+ geom_tile(aes(x = x, y = y, fill = Density)) +
+ coord_equal() +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ scale_fill_viridis(option="C", direction = 1)
+
+
+#### Calculate people within circle for cell centroids, a lot of work to do here
+library(raster)
+
+# subset for testing code
+
+#subset_extent <- ext(4330000, 4380000, 2780000, 2830000)
+
+#hnv_sub <- crop(agg_hnv, subset_extent)
+hnv_sub <- aggregate(hnv, fact=50, fun=mean) # 5km x 5km cells
+
+hnv_sub_gg <- as.data.frame(hnv_sub, xy=TRUE)
+hnv_sub_gg <- rowid_to_column(hnv_sub_gg, var = "ID")
+hnv_to_merge <- hnv_sub_gg %>% dplyr::select(ID, "HNV_share" = "eea_r_3035_100_m_hnv-farmland-ac_2012_v1_r00")
+
+#### 9 Calculate WTP per cell ####
+check_list <- c(unlist(values(hnv_sub)))
+
+agg_hnv_wtp <- -(model1$Estimate[[4]] + 2*model1$Estimate[5]*(check_list*100))/model1$Estimate[6]
+
+values(hnv_sub) <- agg_hnv_wtp/scale_factor
+
+# transform density raster to same resolution #
+dens_sub <- aggregate(density_raster, fact=20, fun=sum) # 5km x 5km cells
+
+grid <- rast(hnv_sub, vals = values(hnv_sub)) # empty raster
+
+grid_points <- as.points(grid, values=TRUE, na.rm=TRUE, na.all=FALSE)
+
+grid_points_gg <- as.data.frame(grid, xy=TRUE)
+
+
+gridPoints <- SpatialPoints(grid_points_gg)
+
+vec_gP <- vect(gridPoints)
+crs(vec_gP) <- crs(dens_sub)
+
+
+buffer_pts <- buffer(vec_gP, buffer_dist)
+
+# Compute people within Buffer
+density_in_buffer <- extract(dens_sub, buffer_pts, fun= sum, na.rm =T, weights = TRUE)
+
+# Check if this is really correct with the ID
+
+grid_points_gg <- rowid_to_column(grid_points_gg, var = "ID")
+
+grid_merge <- left_join(grid_points_gg, density_in_buffer, by="ID")
+
+colnames(grid_merge) <- c("ID", "x", "y", "marg_WTP", "Density")
+
+grid_merge <- left_join(grid_merge, hnv_to_merge, by="ID")
+
+## calculate WTP for 1ha increase, produces infinity values if SQ = 0 ! ##
+#### work on here ### 1% increase is always the same area actually not depended on SQ with the current design
+#need to define this clearly for the DCE #
+# For now assumed that we value 1% increase of HNV on the total land area which is always the same area
+# Aggregated WTP is the Aggregated WTP for a one hectare increase of HNV in the raster cell
+grid_merge <- grid_merge %>% mutate(SQ_HNV = (HNV_share * unique(res(hnv_sub))^2)/10000,
+ agg_WTP = marg_WTP*Density/(((unique(res(hnv_sub))^2)/10000)*0.01))
+
+###
+ggplot(data=grid_merge) +
+ geom_tile(aes(x = x, y = y, fill = agg_WTP)) +
+ coord_equal() +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ scale_fill_viridis(option="C", direction = 1)
+
+ggsave("Figures/agg_wtp_example.png", dpi="print", width = 7, height = 5)
+
+
+# discontinuous map with cuts
+grid_merge$cuts <- cut(grid_merge$agg_WTP, c(-1000, 0, 100, 200, 350, 1000, 3000, 5000, 8600))
+
+
+
+ggplot(data=grid_merge) +
+ geom_tile(aes(x = x, y = y, fill = cuts)) +
+ coord_equal() +
+ annotation_scale(location = "bl", width_hint = 0.1,plot_unit="m", line_width = 0.1) +
+ annotation_north_arrow(location = "br" ,pad_x = unit(0, "in"), pad_y = unit(0, "in"),
+ style = north_arrow_fancy_orienteering) +
+ labs(fill="Aggregated WTP (€/ha)") +
+ xlab("") +
+ ylab("") +
+ scale_fill_viridis(option="C", direction = 1, discrete = T,
+ labels=c("< 0", "0 - 100", "100 - 200", "200 - 350", "350 - 1,000", "1,000 - 3,000", "3,000 - 5,000",
+ " > 5,000"))
+
+ggsave("Figures/agg_wtp_hnv_ex_disc.png", dpi="print", width=6, height = 5)
+
+
+
+
+
diff --git a/Projects/ValuGaps/designsyntax/design syntax.ngs b/Projects/ValuGaps/designsyntax/design syntax.ngs
new file mode 100644
index 0000000..e64339c
--- /dev/null
+++ b/Projects/ValuGaps/designsyntax/design syntax.ngs
@@ -0,0 +1,234 @@
+? Orthogonal Design
+
+
+
+
+design
+;alts = alt1*, alt2*, alt3
+;rows = 72
+;block = 6
+;orth = sim
+
+
+;con
+;model:
+U(alt1) = b1 + b2 * B[0,1,2,3] + b3 * C[0,1,2,3] + b5 * P[40,80,120,160]/
+U(alt2) = b1 + b2 * B + b3 * C + b5 * P
+
+
+
+
+;formatTitle = 'Scenario <scenarionumber> Block <blocknumber>'
+;formatTableDimensions = 3, 5
+;formatTable:
+1,1 = '' /
+1,2 = 'Share of protected areas' /
+1,3 = 'Share of HNV' /
+1,4 = 'Cost per year' /
+1,5 = 'Choice question&:' /
+2,1 = 'alt1' /
+2,2 = '<alt1.b>' /
+2,3 = '<alt1.c>' /
+2,4 = '<alt1.p>' /
+2,5 = '' /
+3,1 = 'alt2' /
+3,2 = '<alt2.b>' /
+3,3 = '<alt2.c>' /
+3,4 = '<alt2.p>' /
+3,5 = ''
+;formatTableStyle:
+1,1 = 'default' /
+1,2 = 'headingattribute' /
+1,3 = 'headingattribute' /
+1,4 = 'headingattribute' /
+1,5 = 'headingattribute' /
+2,1 = 'heading1' /
+2,2 = 'body1' /
+2,3 = 'body1' /
+2,4 = 'body1' /
+2,5 = 'choice1' /
+3,1 = 'heading2' /
+3,2 = 'body2' /
+3,3 = 'body2' /
+3,4 = 'body2' /
+3,5 = 'choice2'
+;formatStyleSheet = Default.css
+;formatAttributes:
+alt1.b(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt1.c(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt1.p(40=#, 80=#, 120=#, 160=#) /
+alt2.b(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt2.c(0=#0% more, 1=#0% more, 2=#0% more, 3=#0% more) /
+alt2.p(40=#, 80=#, 120=#, 160=#)
+$
+
+
+
+
+
+
+
+??? Baysian Efficient Design 1
+
+design
+;alts = alt1*, alt2*, alt3
+;rows = 36
+;block = 3
+;eff = (mnl,d,mean)
+;rep = 1000
+;bdraws = halton(1000)
+;con
+;model:
+U(alt1) = b1[(n,-1.2,0.1)] + b2[(n,0.1,0.02)] * B[0,1] + b3[(n,0.4,0.1)] * C[0,1] + b4[(n,0.3,0.1)] * D[0,1] + b5[(n,0.02,0.001)] * P[35:83:6]/
+U(alt2) = b6[(n,-1.4,0.1)] + b2 * B + b3 * C + b4 * D + b5 * P
+?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)
+$
+
+
+
+?Efficient Design
+
+
+
+design
+;alts = alt1*, alt2*, alt3
+;rows = 36
+;block = 3
+;eff = (mnl,d)
+;rep = 1000
+
+;con
+;model:
+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) =
+
+
+$
+
+
+
+
+
+
+?Bad prior Design
+
+
+design
+;alts = alt1*, alt2*, alt3
+;rows = 18
+;block = 2
+;eff = (mnl,d)
+;rep = 1000
+
+?;con
+;model:
+U(alt1) = b1[-10] + b2[9] * B[0,1] + b3[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
+
+
+;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
new file mode 100644
index 0000000..8e3e7fc
--- /dev/null
+++ b/Projects/ValuGaps/parameters_valugaps.R
@@ -0,0 +1,129 @@
+
+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")
+
+}
--
GitLab