From 3e30039450368785f1e8fa6fd5c2c4058f51051b Mon Sep 17 00:00:00 2001 From: Chris Free Date: Mon, 17 Jan 2022 22:23:09 -0800 Subject: [PATCH] cost search analysis --- ...p1b_reforecast_with_sensitivity_analysis.R | 20 +- ...ivity_analysis_costs_that_change_results.R | 251 ++++++ .../Step2_create_composite_potentials.R | 175 ++++ .../Step3_format_composite_potentials.R | 152 ++++ .../Step4_build_finfish_scenario_results.R | 778 ++++++++++++++++++ .../Step5_build_bivalve_scenario_results.R | 423 ++++++++++ .../Step6_merge_capture_mariculture_results.R | 145 ++++ .../cost_search/figure_cost_search_analysis.R | 72 ++ 8 files changed, 1999 insertions(+), 17 deletions(-) create mode 100644 code/cost_search/Step1c_reforecast_with_sensitivity_analysis_costs_that_change_results.R create mode 100644 code/cost_search/Step2_create_composite_potentials.R create mode 100644 code/cost_search/Step3_format_composite_potentials.R create mode 100644 code/cost_search/Step4_build_finfish_scenario_results.R create mode 100644 code/cost_search/Step5_build_bivalve_scenario_results.R create mode 100644 code/cost_search/Step6_merge_capture_mariculture_results.R create mode 100644 code/cost_search/figure_cost_search_analysis.R diff --git a/code/Step1b_reforecast_with_sensitivity_analysis.R b/code/Step1b_reforecast_with_sensitivity_analysis.R index 065e972..d28292b 100644 --- a/code/Step1b_reforecast_with_sensitivity_analysis.R +++ b/code/Step1b_reforecast_with_sensitivity_analysis.R @@ -43,9 +43,9 @@ data <- data %>% rcp <- "rcp26" species <- data[1,] price_scalar <- 0.7 -cost_scalar <- 1.2 +cost_scalar <- 1.3 periods <- c("2021-2030", "2051-2060", "2091-2100") -reforecast <- function(species, rcp, periods, price_scalar=1, cost_scalar=1, outdir=F, plot=T){ +reforecast <- function(species, rcp, periods, price_scalar=0.7, cost_scalar=1.3, outdir=F, plot=T){ # 1. Read data ############################ @@ -219,23 +219,8 @@ for(i in 1:nrow(data)){ species <- data[i,] periods <- c("2021-2030", "2051-2060", "2091-2100") output <- reforecast(species=species, periods=periods, rcp="rcp26", outdir=outputdir, plot=T) -} - -for(i in 1:nrow(data)){ - species <- data[i,] - periods <- c("2021-2030", "2051-2060", "2091-2100") output <- reforecast(species=species, periods=periods, rcp="rcp45", outdir=outputdir, plot=F) -} - -for(i in 1:nrow(data)){ - species <- data[i,] - periods <- c("2021-2030", "2051-2060", "2091-2100") output <- reforecast(species=species, periods=periods, rcp="rcp60", outdir=outputdir, plot=F) -} - -for(i in 1:nrow(data)){ - species <- data[i,] - periods <- c("2021-2030", "2051-2060", "2091-2100") output <- reforecast(species=species, periods=periods, rcp="rcp85", outdir=outputdir, plot=F) } @@ -245,3 +230,4 @@ for(i in 1:nrow(data)){ + diff --git a/code/cost_search/Step1c_reforecast_with_sensitivity_analysis_costs_that_change_results.R b/code/cost_search/Step1c_reforecast_with_sensitivity_analysis_costs_that_change_results.R new file mode 100644 index 0000000..4360b30 --- /dev/null +++ b/code/cost_search/Step1c_reforecast_with_sensitivity_analysis_costs_that_change_results.R @@ -0,0 +1,251 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(raster) +library(ggplot2) +library(tidyverse) + +# Directories +codedir <- "code" +sppdir <- "data/species/data" +indir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/raw" +outdir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search" +plotdir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_plots" + +# Read aquacast function +source(file.path(codedir, "calc_costs_v2.R")) + +# Read species data +load(file.path(sppdir, "aquaculture_species_key_20cages.Rdata")) + +# Brackish ISSCAAPs +brackish_isscaaps <- c("Freshwater molluscs", "Miscellaneous diadromous fishes", + "Miscellaneous freshwater fish", "River eels", "Shads", + "Sturgeons, paddlefishes", "Tilapias and other cichlids") + +# Not lonline bivalves +bad_bivalve_isscaaps <- c("Clams, cockles, arkshells", "Pearls, mother-of-pearl, shells", "Scallops, pectens") + +# Reduce data +data <- data %>% + filter(!isscaap %in% c(brackish_isscaaps, bad_bivalve_isscaaps)) + + + +# Build function +################################################################################ + +# Function +rcp <- "rcp26" +species <- data[1,] +cost_scalar <- 1.2 +periods <- c("2021-2030", "2051-2060", "2091-2100") +reforecast <- function(species, rcp, periods, price_scalar=1, cost_scalar=1, outdir=F, plot=T){ + + # 1. Read data + ############################ + + # Param + spp <- species$species + type <- ifelse(species$class=="Bivalvia", "Bivalve", "Finfish") + print(spp) + + # Read suitable cells + print("... reading and transforming habitat suitability") + filename_base <- paste(toupper(rcp), gsub(" ", "_", spp), sep="_") + filename_in <- paste0(filename_base, ".Rds") + vcells_orig <- readRDS(file.path(indir, filename_in)) + periods_avail <- sort(unique(vcells_orig$period)) + + # Convert to raster + vcells <- vcells_orig %>% + # Simplify + select(period, x, y, viable) %>% + # Spread to XY-ZZZ + spread(key="period", value="viable") %>% + # Convert to raster + rasterFromXYZ(crs = crs(ras_temp)) %>% + # Project raster + projectRaster(to=ras_temp) + + # 2. Extract data + ############################ + + # Growth and harvest parameters + linf_cm <- species$linf_cm + k <- species$k + harvest_cm <- species$harvest_cm + harvest_g <- species$harvest_g + harvest_yr <- species$harvest_yr + harvest_kg_m3 <- species$harvest_kg_m3 + harvest_cm_ft <- species$harvest_cm_ft + nstocked <- species$nstocked + nstocked_adj <- species$nstocked_adj + a <- species$a + b <- species$b + fcr <- species$fcr + price_usd_mt <- species$price_usd_mt_isscaap * price_scalar + + + # 3. Farm design + ############################ + + # Finfish farm design + if(type=="Finfish"){ + # Finfish farm design + farm_design <- tibble(type="finfish", + area_sqkm=1, + ncages=20, + cage_vol_m3=9000) %>% + # Calculate number stocked + mutate(tot_m3=ncages * cage_vol_m3, + tot_kg=tot_m3*harvest_kg_m3, + nstocked=tot_kg*1000/harvest_g, + nstocked_adj=nstocked) + } + + # Bivalve farm design + if(type=="Bivalve"){ + farm_design <- tibble(type="bivalve", + area_sqkm=1, + nlines=species$lines_n, + line_rope_ft=2109, + harvest_cm_ft=harvest_cm_ft) %>% + mutate(nstocked = nlines * line_rope_ft * (harvest_cm_ft / harvest_cm), + nstocked_adj=nstocked_adj) + } + + # Check number of stocked + if(all.equal(nstocked, farm_design$nstocked)!=T){ + stop("Number of stocked individuals calculated doesn't match value in key.") + } + + # Calculate yield per year (kg/yr) for 1 sqkm farm + farm_kg_yr <- harvest_g/1000 * farm_design$nstocked / harvest_yr + farm_mt_yr <- farm_kg_yr / 1000 + + # Check farm production + if(all.equal(species$prod_mt_yr, farm_mt_yr)!=T){ + stop("Farm production doesn't match value in key.") + } + + # 4. Production + ############################ + + # Calculate annual harvest per cell + print("... mapping production potential") + cell_sqkm <- 100 + cell_nfarms <- cell_sqkm / farm_design$area_sqkm + cell_prod_mt_yr <- farm_mt_yr * cell_nfarms + + # Calculate annual revenue per cell + cell_revenue_usd_yr <- cell_prod_mt_yr * price_usd_mt + + # Check cell revenue + if(all.equal(species$revenue_usd_yr*cell_nfarms*price_scalar, cell_revenue_usd_yr)!=T){ + stop("Cell revenue doesn't match value in key.") + } + + + # 5. Calculate costs + #################################### + + # Calculate costs and profits + print("... calculating costs and profits") + cell_cost_usd_yr <- calc_costs_v2(farm_design, cell_prod_mt_yr, fcr, vcells, harvest_yr) * cost_scalar + + # 6. Calculate profits + #################################### + + # Build final data frame + data_df <- cell_cost_usd_yr %>% + # Convert to dataframe + as.data.frame(xy=T) %>% + setNames(c("x", "y", periods_avail)) %>% + gather(key="period", value="cost_usd_yr", 3:ncol(.)) %>% + # Reduce to cells with costs + filter(!is.na(cost_usd_yr)) %>% + # Add variables + mutate(viable=1, + prod_mt_yr=cell_prod_mt_yr, + revenue_usd_yr=cell_revenue_usd_yr, + profits_usd_yr=revenue_usd_yr - cost_usd_yr) %>% + # Arrange + select(period, x, y, viable, prod_mt_yr, revenue_usd_yr, cost_usd_yr, profits_usd_yr) + + + # 5. Export and plot + #################################### + + # If exporting + outfile_basename <- paste(toupper(rcp), gsub(" ", "_", spp), "cost_search", sep="_", cost_scalar) + if(outdir!=F){ + saveRDS(data_df, file.path(outdir, paste0(outfile_basename, ".Rds"))) + } + + # If plotting + if(plot==T){ + + # Build stats data frame + stats <- data_df %>% + filter(profits_usd_yr>0) %>% + group_by(period) %>% + summarize(area_sqkm_mil=sum(viable)*100/1e6, + prod_mt_yr_mil=sum(prod_mt_yr)/1e6, + profits_usd_yr_tril=sum(profits_usd_yr)/1e12) %>% + ungroup() + + # Plot + g <- ggplot(stats, aes(x=period, y=area_sqkm_mil)) + + geom_bar(stat="identity") + + labs(x="", y="Profitable area (millions of sqkm)", title=species$comm_name) + + theme_bw() + + theme( axis.text.y = element_text(angle = 90, hjust = 0.5) ) + print(g) + + } + +} + + + +# Reforecast using function +################################################################################ + +# Cost scalars +cost_scalars <- c(2.1, 2.7) + +# Build key +spp_cost_key <-purrr::map_df(cost_scalars, function(x){ + + # Add cost scalar + data1 <- data %>% + mutate(cost_scalar=x) + +}) + + +# Loop through species and cost scalars +for(i in 1:nrow(spp_cost_key)){ + species <- spp_cost_key[i,] + cost_scalar <- spp_cost_key$cost_scalar[i] + periods <- c("2021-2030", "2051-2060", "2091-2100") + output <- reforecast(species=species, periods=periods, rcp="rcp26", cost_scalar=cost_scalar, outdir=outdir, plot=T) + output <- reforecast(species=species, periods=periods, rcp="rcp45", cost_scalar=cost_scalar, outdir=outdir, plot=F) + output <- reforecast(species=species, periods=periods, rcp="rcp60", cost_scalar=cost_scalar, outdir=outdir, plot=F) + output <- reforecast(species=species, periods=periods, rcp="rcp85", cost_scalar=cost_scalar, outdir=outdir, plot=F) +} + + + + + + + + + diff --git a/code/cost_search/Step2_create_composite_potentials.R b/code/cost_search/Step2_create_composite_potentials.R new file mode 100644 index 0000000..9a1caf0 --- /dev/null +++ b/code/cost_search/Step2_create_composite_potentials.R @@ -0,0 +1,175 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(raster) +library(ggplot2) +library(tidyverse) + +# Directories +sppdir <- "data/species/data" +inputdir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search" +outputdir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_processed" + +# Read species data +load(file.path(sppdir, "aquaculture_species_key.Rdata")) +spp_key_orig_full <- data +rm(data_full) + +# Brackish ISSCAAPs +brackish_isscaaps <- c("Freshwater molluscs", "Miscellaneous diadromous fishes", + "Miscellaneous freshwater fish", "River eels", "Shads", + "Sturgeons, paddlefishes", "Tilapias and other cichlids") + +# Not lonline bivalves +bad_bivalve_isscaaps <- c("Clams, cockles, arkshells", "Pearls, mother-of-pearl, shells", "Scallops, pectens") + +# Remove species from species key +spp_key_orig <- spp_key_orig_full %>% + filter(!isscaap %in% c(brackish_isscaaps, bad_bivalve_isscaaps)) + + +# Setup +################################################################################ + +# Function to merge results +rcp <- "rcp85"; type <- "bivalve"; outdir <- outputdir; suffix <- "cost_search_1.3" +merge_results <- function(rcp, type, outdir, suffix=""){ + + # Species key + type_do <- ifelse(type=="finfish", "Actinopterygii", "Bivalvia") + if(suffix==""){ + spp_key <- spp_key_orig %>% + filter(class==type_do) %>% + select(species, comm_name) %>% + mutate(file_name=paste0(toupper(rcp), "_", gsub(" ", "_", species), ".Rds")) + }else{ + spp_key <- spp_key_orig %>% + filter(class==type_do) %>% + select(species, comm_name) %>% + mutate(file_name=paste0(toupper(rcp), "_", gsub(" ", "_", species), "_", suffix, ".Rds")) + } + + # Identify files to merge + files_do <- spp_key$file_name + + # Chunk files into groups of 20 + # I have to do this b/c I was using up all memory on finfish and crashing R + # x <- 1 + files_do_chunks <- split(files_do, ceiling(seq_along(files_do)/20)) + data <- purrr::map_df(1:length(files_do_chunks), function(x){ + + # Chunk to do + files_do_chunk <- files_do_chunks[[x]] + + # Loop through files and merge + # y <- files_do_chunk[1] + data_chunk <- purrr:::map_df(files_do_chunk, function(y){ + + # Read one file + # file_do <- files_do[1] + file_do <- y + if(suffix==""){ + spp_do <- file_do %>% gsub(".Rds", "", .) %>% gsub(paste0(toupper(rcp), "_"), "", .) %>% gsub("_", " ", .) + }else{ + spp_do <- file_do %>% gsub(suffix, "", .) %>% gsub(".Rds", "", .) %>% gsub(paste0(toupper(rcp), "_"), "", .) %>% gsub("_", " ", .) %>% stringr::str_trim() + } + sdata_orig <- readRDS(file.path(inputdir, file_do)) + sdata <- sdata_orig %>% + filter(profits_usd_yr>0) %>% + mutate(species=spp_do) %>% + select(species, everything()) + + }) + + # Identify most profitable species for a cell in this chunk + results_chunk <- data_chunk %>% + group_by(period, x, y) %>% + arrange(period, x, y, desc(profits_usd_yr)) %>% + slice(1) %>% + ungroup() + + }) + + # Identify most profitable species for a cell across all chunks + results <- data %>% + group_by(period, x, y) %>% + arrange(period, x, y, desc(profits_usd_yr)) %>% + slice(1) %>% + ungroup() + + # Calculate year stats + ystats <- results %>% + group_by(period) %>% + summarise(ncells=n(), + area_sqkm=n()*100/1e6, + prod_mt=sum(prod_mt_yr)/1e9, + profits_usd=sum(profits_usd_yr)/1e12) %>% + ungroup() + + # Plot production trends + g1 <- ggplot(ystats, aes(x=period, y=prod_mt)) + + geom_bar(stat="identity") + + labs(x="", y="Production (billions of mt)") + + theme_bw() + print(g1) + + # Plot 2100 production + period_last <- unique(ystats$period)[length(unique(ystats$period))] + title_text <- paste(toupper(rcp), str_to_title(type), period_last, "production") + g2 <- ggplot(filter(results, period==period_last), + aes(x=x, y=y, fill=prod_mt_yr/1e6)) + + geom_tile() + + labs(x="", y="", title=title_text) + + scale_fill_gradientn(name="Production (millions mt)", + colors=rev(RColorBrewer::brewer.pal(9, "RdBu"))) + + theme_bw() + print(g2) + + # Export data + if(suffix==""){ + outfile <- paste(toupper(rcp), str_to_title(type), "rational.Rds", sep="_") + }else{ + outfile <- paste0(toupper(rcp), "_", str_to_title(type), "_rational_", suffix, ".Rds") + } + saveRDS(results, file.path(outdir, outfile)) + +} + + +# Setup +################################################################################ + +# Run +cost_scalars <- c(2.1, 2.3, 2.5, 2.7) + +# Loop through scalars +for(i in 1:length(cost_scalars)){ + + # Build suffix + suffix <- paste0("cost_search_", cost_scalars[i]) + + # Bivalves + merge_results(rcp="RCP26", type="bivalve", outdir=outputdir, suffix=suffix) + merge_results(rcp="RCP45", type="bivalve", outdir=outputdir, suffix=suffix) + merge_results(rcp="RCP60", type="bivalve", outdir=outputdir, suffix=suffix) + merge_results(rcp="RCP85", type="bivalve", outdir=outputdir, suffix=suffix) + + # Finfish - 11 minutes each + merge_results(rcp="RCP26", type="finfish", outdir=outputdir, suffix=suffix) + merge_results(rcp="RCP45", type="finfish", outdir=outputdir, suffix=suffix) + merge_results(rcp="RCP60", type="finfish", outdir=outputdir, suffix=suffix) + merge_results(rcp="RCP85", type="finfish", outdir=outputdir, suffix=suffix) + +} + + + + + + + diff --git a/code/cost_search/Step3_format_composite_potentials.R b/code/cost_search/Step3_format_composite_potentials.R new file mode 100644 index 0000000..3f6431e --- /dev/null +++ b/code/cost_search/Step3_format_composite_potentials.R @@ -0,0 +1,152 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(raster) +library(ggplot2) +library(tidyverse) +library(countrycode) +library(fuzzyjoin) + +# Directories +eezdir <- "data/eezs" +sppdir <- "data/species/data" +maskdir <- "data/constraints/masks" +datadir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_processed" + +# Read EEZ data +eezs <- raster(file.path(eezdir, "eezs_v10_raster_10km.tif")) +eez_key_orig <- read.csv(file.path(eezdir, "eezs_v10_key.csv"), as.is=T) + +# Read, format, and check mask +mask <- raster(file.path(maskdir, "mask_10km.grd")) +mask[mask==0] <- NA +plot(mask) + +# Convert mask to data frame +mask_df <- mask %>% + # Convert mask to dataframe + as.data.frame(xy=T) %>% + setNames(c("x", "y", "use")) %>% + # Reduce to viable cells + filter(!is.na(use)) %>% + # Fix coordinates (needed for successful merge) + mutate(x=as.integer(x), + y=as.integer(y)) + +# Format EEZ data +################################################################################ + +# Continent keys +south_ind <- c("Amsterdam and Saint Paul Islands", "Heard and McDonald Islands", + "Kerguélen", "Cocos Islands", "Chagos Archipelago", "Crozet Islands") +south_atl <- c("Ascension", "Tristan da Cunha", "South Georgia and the South Sandwich Islands", "Bouvet") +eq_pacific <- c("Howland and Baker islands", "Palmyra Atoll", "Wake Island", "Jarvis Island", "Johnston Atoll") +moz_channel <- c("Europa Island", "Juan de Nova Island", "Bassas da India") + +# Build EEZ key +eez_key <- eez_key_orig %>% + # Add continent + mutate(ter1_continent=countrycode(sourcevar = ter1_iso, origin="iso3c", destination = "continent"), + sov1_continent=countrycode(sourcevar = sov1_iso, origin="iso3c", destination = "continent")) %>% + # Fill missing continent typ + mutate(ter1_continent=ifelse(ter1_name %in% south_ind, "Indian Ocean", ter1_continent), + ter1_continent=ifelse(ter1_name %in% south_atl, "S Atlantic Ocean", ter1_continent), + ter1_continent=ifelse(ter1_name %in% eq_pacific, "Oceania", ter1_continent), + ter1_continent=ifelse(ter1_name %in% moz_channel, "Africa", ter1_continent), + ter1_continent=ifelse(ter1_name == "Clipperton Island", "Americas", ter1_continent), + ter1_continent=ifelse(ter1_name == "Antarctica", "Antarctica", ter1_continent), + sov1_continent=ifelse(is.na(sov1_continent), ter1_continent, sov1_continent)) + +# Convert EEZ raster to dataframe +eezs_df <- as.data.frame(eezs, xy=T) %>% + setNames(c("x", "y", "eez_code")) %>% + left_join(select(eez_key, eez_code, eez_name, eez_type, + ter1_name, ter1_iso, sov1_name, sov1_iso, + ter1_continent, sov1_continent)) %>% + # Convert coords to integers to ease left_join below + mutate(x=as.integer(x), + y=as.integer(y)) + + +# Add EEZ data +################################################################################ + +# Function +# rcp <- "rcp26"; type <- "finfish"; suffix="new_costs" +format_data <- function(rcp, type, suffix=""){ + + # Read data + if(suffix==""){ + infile <- paste(toupper(rcp), str_to_title(type), "rational.Rds", sep="_") + outfile <- paste(toupper(rcp), str_to_title(type), "rational_use.Rds", sep="_") + }else{ + infile <- paste0(toupper(rcp), "_", str_to_title(type), "_rational_", suffix, ".Rds") + outfile <- paste0(toupper(rcp), "_", str_to_title(type), "_rational_use_", suffix, ".Rds") + } + data_orig <- readRDS(file.path(datadir, infile)) + + # Format data + data <- data_orig %>% + # Format xy for join + mutate(x=as.integer(x), + y=as.integer(y)) %>% + # Add mask and reduce to viable areas + left_join(mask_df, by=c("x", "y")) %>% + filter(use==1) %>% + select(-use) %>% + # Add EEZ information + left_join(eezs_df, by=c("x", "y")) %>% + # Remove cells with out SOV1 info (pertain to an overlapping claim) + filter(!is.na(sov1_name)) + + # Inspect + freeR::complete(data) + + # Summarize plot + stats <- data %>% + group_by(period, ter1_continent) %>% + summarize(prod_bt_yr=sum(prod_mt_yr)/1e9) + g <- ggplot(stats, aes(x=period, y=prod_bt_yr, fill=ter1_continent)) + + geom_bar(stat="identity") + + labs(x="Period", y="Production potential (billions mt)") + + theme_bw() + print(g) + + # Export data + saveRDS(data, file.path(datadir, outfile)) + +} + +# Run +cost_scalars <- c(2.1, 2.3, 2.5, 2.7) + +# Loop through scalars +for(i in 1:length(cost_scalars)){ + + # Build suffix + suffix <- paste0("cost_search_", cost_scalars[i]) + + # Finfish + format_data(rcp="rcp26", type="finfish", suffix=suffix) + format_data(rcp="rcp45", type="finfish", suffix=suffix) + format_data(rcp="rcp60", type="finfish", suffix=suffix) + format_data(rcp="rcp85", type="finfish", suffix=suffix) + + # Bivalves + format_data(rcp="rcp26", type="bivalve", suffix=suffix) + format_data(rcp="rcp45", type="bivalve", suffix=suffix) + format_data(rcp="rcp60", type="bivalve", suffix=suffix) + format_data(rcp="rcp85", type="bivalve", suffix=suffix) + +} + + + + + + diff --git a/code/cost_search/Step4_build_finfish_scenario_results.R b/code/cost_search/Step4_build_finfish_scenario_results.R new file mode 100644 index 0000000..1fdfdeb --- /dev/null +++ b/code/cost_search/Step4_build_finfish_scenario_results.R @@ -0,0 +1,778 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(raster) +library(ggplot2) +library(tidyverse) +library(countrycode) +library(rnaturalearth) + +# Directories +datadir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/processed" +outputdir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_processed" +plotdir <- "figures" +tabledir <- "tables" + +# Read species data +load("data/species/data/aquaculture_species_key.Rdata") +spp_key <- data +rm(data_full, data) + +# Read feed info +fifos <- read.csv("tables/TableS6_fifo_fcr_fmfo_stats.csv", as.is=T) +ffdata <- read.csv("data/feed_params/processed/forage_fish_availability.csv", as.is=T) +maq_prod <- read.csv("data/feed_params/processed/FAO_2013_2017_maq_prod_averages_by_country.csv", as.is=T) + +# World +world <- rnaturalearth::ne_countries(scale="large", type = "countries", returnclass = "sf") + +# Read EEZ key +eezdir <- "data/eezs" +eez_key <- read.csv(file.path(eezdir, "eezs_v10_key.csv"), as.is=T) + +# Read population data +popdir <- "data/pop_growth/data" +pop <- readRDS(file.path(popdir, "WB_UN_1960_2100_human_population_by_country.Rds")) + +# Read fisheries impacts +wc <- readRDS("data/capture_projections/data/Free_etal_2020_national_projections_with_pop_data.Rds") + + +# Historical data +################################### + +# Read historical data +hist_orig <- readRDS(file.path(datadir, "FAO_1950_2018_wc_aq_seafood_per_capita_national.Rds")) + +# Build today's values +hist <- hist_orig %>% + # Recent year + filter(year==max(year)) %>% + # Sum across sectors + group_by(country, iso3, npeople) %>% + summarize(prod_mt=sum(prod_mt), + meat_mt=sum(meat_mt), + meat_kg_person_2017=sum(meat_kg_person)) %>% + ungroup() + + +# Function to develop feed-constrained finfish mariculture +################################################################################ + +# For testing +rcp <- "RCP 2.6" +mgmt_scenario <- "Reformed fisheries management" +feed_scenario <- "Reformed feed use" +dev_scenario <- "Optimum" +suffix <- "cost_search_3" + +# Function to develop feed-constrained finfish mariculture +expand_mariculture <- function(rcp, mgmt_scenario, feed_scenario, dev_scenario, suffix=""){ + + # Steps + # 1. Determine feed availability + # 2. Determine feed demand + # 3. Develop based on distribution pattern + + # Setup + ############################################################### + + # Parameters + rcp_do <- rcp + mgmt_scenario_do <- mgmt_scenario + feed_scenario_do <- feed_scenario + dev_scenario_do <- dev_scenario + + # Demand limit + prod_mt_cap <- 103 * 0.27 / 0.87 * 1e6 + + # Read finfish forecasts + rcp_do_short <- gsub("\\.| ", "", rcp_do) + if(suffix==""){ + infile <- paste0( rcp_do_short, "_Finfish_rational_use.Rds") + }else{ + infile <- paste0( rcp_do_short, "_Finfish_rational_use_", suffix, ".Rds") + } + data_orig <- readRDS(file.path(outputdir, infile)) + + # Build forecasts period key + periods <- data_orig %>% pull(period) %>% unique() + period_key <- purrr::map_df(periods, function(x){ + yr1 <- substr(x,1,4) %>% as.numeric() + yr2 <- substr(x,6,9) %>% as.numeric() + yrs <- yr1:yr2 + yr_df <- tibble(period=x, + year=yrs) + }) + + # Determine feed availability + ############################################################### + + # Determine forage fish availability given RCP, fisheries management, and feed use scenario + ff_avail <- ffdata %>% + # Reduce to RCP and fisheries management scenario + filter(rcp==rcp_do & mgmt_scenario==mgmt_scenario_do & feed_scenario==feed_scenario_do) %>% + # Add forecast period + left_join(period_key, by="year") %>% + # Simplify and rename columns + select(rcp, period, year, catch_ff_mt_maq) %>% + rename(ff_avail_mt_yr=catch_ff_mt_maq) %>% + # Remove years not in period + filter(!is.na(period)) %>% + # Calculate mean over forecast periods + group_by(rcp, period) %>% + summarize(ff_avail_mt_yr=mean(ff_avail_mt_yr)) %>% + ungroup() + + # Determine feed demands + ############################################################### + + # Determine which FIFO ratios to use based on feed scenario + fifos_use <- fifos + if(feed_scenario_do=="Reformed feed use"){ + fifos_use$fifo_use <- fifos_use$fifo2050 + }else{ + fifos_use$fifo_use <- fifos_use$fifo2030 + } + + # Determine feed demands of each profitable cell (based on feed use scenario) + data1 <- data_orig %>% + # Add feed group and FIFO ratios + left_join(select(spp_key, species, feed_group), by="species") %>% + left_join(select(fifos_use, group, fifo_use), by=c("feed_group"="group")) %>% + # Add forage fish availability + left_join(ff_avail, by=c("period")) %>% + # Calculate forage fish demand + mutate(ff_demand_mt_yr=prod_mt_yr*fifo_use) %>% + # Make sure Alaska and Hawaii have US designation + mutate(ter1_name_use=recode(ter1_name, + "Hawaii"="United States", + "Alaska"="United States", + "Micronesia"="Federated States of Micronesia", + "Comores"="Comoros"), + ter1_iso_use=countrycode(ter1_name_use, "country.name", "iso3c")) %>% + # Use the sovereign nation as the jurisdictional unit for the remaining territories without ISOs + mutate(ter1_name_use=ifelse(is.na(ter1_iso_use), sov1_name, ter1_name_use)) %>% + # Run ISO3 lookup one last time and now it should be perfect + mutate(ter1_iso_use=countrycode(ter1_name_use, "country.name", "iso3c")) + + + # Develop mariculture based on development scenario + ############################################################### + + # Rational development + if(dev_scenario_do=="Rational"){ + + # Develop profitable cells first + data2 <- data1 %>% + group_by(period) %>% + arrange(period, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish + mutate(cum_ff_demand_mt_yr=cumsum(ff_demand_mt_yr), + cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_ff_demand_mt_yr <= ff_avail_mt_yr & cum_prod_mt_yr <= prod_mt_cap, "yes", "no")) %>% + filter(developed=="yes") + + } + + # Current development + if(dev_scenario_do=="Current"){ + + # Format current producers + # Based on FAO 2013-2017 averages + props <- maq_prod %>% + # Reduce to finfish aquaculture totals + filter(major_group=="Pisces") %>% + # Format country information + mutate(country_use=countrycode(country_use, "country.name", "country.name"), + iso3_use=countrycode(country_use, "country.name", "iso3c")) %>% + # Remove TERRITORIES not forecast to have to viable mariculutre + filter(iso3_use %in% unique(data1$ter1_iso_use)) %>% + # Calculate proportion of feed dedicated to this TERRITORY + mutate(feed_prop=prod_mt/sum(prod_mt)) # note if a country isn't present in a period some feed will go unused + + # Develop profitable cells first, in countries with production + data2 <- data1 %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3_use, feed_prop), by=c("ter1_iso_use"="iso3_use")) %>% + # Calculate forage fish allocation + mutate(feed_prop=ifelse(is.na(feed_prop), 0, feed_prop), + feed_allocated_mt_yr=ff_avail_mt_yr*feed_prop) %>% + # Remove places not allocated any forage fish + filter(feed_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_ff_demand_mt_yr=cumsum(ff_demand_mt_yr), + developed=ifelse(cum_ff_demand_mt_yr<=feed_allocated_mt_yr, "yes", "no")) %>% + filter(developed=="yes") + + # If feed isn't limiting (i.e., demand is exceeded), distribute as proportion of total global demand + # NOTE: THIS CODE IS IMPERFECT. IT WILL FAIL IF ONE PERIOD IS BELOW AND ANOTHER IS ABOVE. _ WROT E A STOP TO FLAG IF HAPPENS + prod_my_yr_period <- data2 %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) + cap_exceeded_yn <- sum(prod_my_yr_period$prod_mt_yr > prod_mt_cap) + if(cap_exceeded_yn%in%c(1,2)){stop("WRONG")} + if(cap_exceeded_yn==3){ + + # Pass #1 + pass1 <- data1 %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3_use, feed_prop), by=c("ter1_iso_use"="iso3_use")) %>% + # Calculate production allocation + rename(prod_prop=feed_prop) %>% + mutate(prod_prop=ifelse(is.na(prod_prop), 0, prod_prop), + prod_allocated_mt_yr=prod_mt_cap*prod_prop) %>% + # Remove places not allocated any forage fish + filter(prod_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_allocated_mt_yr, "yes", "no")) + + # Cells developed in pass one + pass1_use <- pass1 %>% + filter(developed=="yes") + + # Undeveloped production from pass 1 + pass1_leftover <- pass1_use %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) %>% + ungroup() %>% + mutate(prod_mt_yr_fill=prod_mt_cap-prod_mt_yr) + + # Pass #2: Develop current countries to fill remaining demand from the pass 1 + pass2_use <- pass1 %>% + # Reduce to undeveloped cells + filter(developed=="no") %>% + # Add remaining production available per period + left_join(pass1_leftover %>% select(period, prod_mt_yr_fill), by="period") %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_yr_fill, "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge the two passes + data2 <- bind_rows(pass1_use, pass2_use) + + + } + + } + + # Proportional development + if(dev_scenario_do=="Proportional"){ + + # 2050 population + pop55 <- pop %>% + filter(year==2100) + + # ISO key for production forecasts + props <- data1 %>% + select(sov1_iso, sov1_name, ter1_iso_use, ter1_name_use) %>% + unique() %>% + arrange(ter1_iso_use) %>% + # Add 2055 population estimates + left_join(pop55 %>% select(iso3, pop_size_50perc), by=c("ter1_iso_use"="iso3")) %>% + # Only consider areas with population projection estimates + filter(!is.na(pop_size_50perc)) %>% + # Calculate proportional allocation + mutate(feed_prop=pop_size_50perc/sum(pop_size_50perc)) %>% + # Simplify + select(ter1_iso_use, feed_prop) + + # Plot allocations + #hist(props$feed_prop, breaks=seq(0, 0.25, 0.005)) + #abline(v=1/nrow(props)) + + # Develop profitable cells first, in countries with production + data2 <- data1 %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(props, by=c("ter1_iso_use")) %>% + # Calculate forage fish allocation + mutate(feed_prop=ifelse(is.na(feed_prop), 0, feed_prop), + feed_allocated_mt_yr=ff_avail_mt_yr*feed_prop) %>% + # Remove places not allocated any forage fish + filter(feed_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_ff_demand_mt_yr=cumsum(ff_demand_mt_yr), + developed=ifelse(cum_ff_demand_mt_yr<=feed_allocated_mt_yr, "yes", "no")) %>% + filter(developed=="yes") + + # If feed isn't limiting (i.e., demand is exceeded), distribute as proportion of total global demand + prod_my_yr_period <- data2 %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) + periods_cap_exceeded <- periods[prod_my_yr_period$prod_mt_yr > prod_mt_cap] + periods_cap_not_exceeded <- periods[!periods %in% periods_cap_exceeded] + + # If any periods are above demand + if(length(periods_cap_exceeded)>0){ + + # Data where not exceeded + if(length(periods_cap_not_exceeded)==0){ + # Create empty dataframe with right columns + data_not_exceeded <- data2 %>% slice(0) + }else{ + data_not_exceeded <- data1 %>% + filter(period %in% periods_cap_not_exceeded) %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(props, by=c("ter1_iso_use")) %>% + # Calculate forage fish allocation + mutate(feed_prop=ifelse(is.na(feed_prop), 0, feed_prop), + feed_allocated_mt_yr=ff_avail_mt_yr*feed_prop) %>% + # Remove places not allocated any forage fish + filter(feed_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_ff_demand_mt_yr=cumsum(ff_demand_mt_yr), + developed=ifelse(cum_ff_demand_mt_yr<=feed_allocated_mt_yr, "yes", "no")) %>% + filter(developed=="yes") + } + + + # Data where exceeded + data_exceeded <- data1 %>% + filter(period %in% periods_cap_exceeded) + + # Pass #1 + pass1 <- data_exceeded %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(props, by="ter1_iso_use") %>% + # Calculate production allocation + rename(prod_prop=feed_prop) %>% + mutate(prod_prop=ifelse(is.na(prod_prop), 0, prod_prop), + prod_allocated_mt_yr=prod_mt_cap*prod_prop) %>% + # Remove places not allocated any forage fish + filter(prod_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_allocated_mt_yr, "yes", "no")) + + # Cells developed in pass one + pass1_use <- pass1 %>% + filter(developed=="yes") + + # Undeveloped production from pass 1 + pass1_leftover <- pass1_use %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) %>% + ungroup() %>% + mutate(prod_mt_yr_fill=prod_mt_cap-prod_mt_yr) + + # Pass #2: Develop current countries to fill remaining demand from the pass 1 + pass2_use <- pass1 %>% + # Reduce to undeveloped cells + filter(developed=="no") %>% + # Add remaining production available per period + left_join(pass1_leftover %>% select(period, prod_mt_yr_fill), by="period") %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_yr_fill, "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge the two passes + data2 <- bind_rows(data_not_exceeded, pass1_use, pass2_use) + + } + + } + + # Need-based development: ATTEMPT 1 + if(dev_scenario_do=="Need-based"){ + + # Build proportions + props <- wc %>% + # Reduce to climate scenario + filter(rcp==rcp_do & scenario=="Full Adaptation") %>% + # Calculate amount of seafood needed to break even + group_by(country, iso3) %>% + summarize(meat_mt_need = meat_kg_person[year==2012]/1000*npeople[year==2100], + meat_mt_avail = meat_mt[year==2100], + meat_mt_miss = meat_mt_need - meat_mt_avail) %>% + ungroup() %>% + # Remove countries gaining seafood/capita + filter(meat_mt_miss>0 * !is.na(meat_mt_miss)) %>% + # Remove countries without aquaculture potential + filter(iso3 %in% unique(data_orig$ter1_iso)) %>% + # Calculate proportion of feed dedicated to this TERRITORY + mutate(feed_prop=meat_mt_miss/sum(meat_mt_miss)) + + # Develop profitable cells first, in countries with production + data2 <- data1 %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3, feed_prop), by=c("ter1_iso_use"="iso3")) %>% + # Calculate forage fish allocation + mutate(feed_prop=ifelse(is.na(feed_prop), 0, feed_prop), + feed_allocated_mt_yr=ff_avail_mt_yr*feed_prop) %>% + # Remove places not allocated any forage fish + filter(feed_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_ff_demand_mt_yr=cumsum(ff_demand_mt_yr), + developed=ifelse(cum_ff_demand_mt_yr<=feed_allocated_mt_yr, "yes", "no")) %>% + filter(developed=="yes") + + # If feed isn't limiting (i.e., demand is exceeded), distribute as proportion of total global demand + # NOTE: THIS CODE IS IMPERFECT. IT WILL FAIL IF ONE PERIOD IS BELOW AND ANOTHER IS ABOVE. _ WROT E A STOP TO FLAG IF HAPPENS + prod_my_yr_period <- data2 %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) + cap_exceeded_yn <- sum(prod_my_yr_period$prod_mt_yr > prod_mt_cap) + if(cap_exceeded_yn%in%c(1,2)){stop("WRONG")} + if(cap_exceeded_yn==3){ + + # Pass #1 + pass1 <- data1 %>% + # Add percent of forage fish allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3, feed_prop), by=c("ter1_iso_use"="iso3")) %>% + # Calculate production allocation + rename(prod_prop=feed_prop) %>% + mutate(prod_prop=ifelse(is.na(prod_prop), 0, prod_prop), + prod_allocated_mt_yr=prod_mt_cap*prod_prop) %>% + # Remove places not allocated any forage fish + filter(prod_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name_use) %>% + arrange(period, ter1_name_use, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative forage fish and see if it exceeds ALLOCATED AMOUNT + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_allocated_mt_yr, "yes", "no")) + + # Cells developed in pass one + pass1_use <- pass1 %>% + filter(developed=="yes") + + # Undeveloped production from pass 1 + pass1_leftover <- pass1_use %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) %>% + ungroup() %>% + mutate(prod_mt_yr_fill=prod_mt_cap-prod_mt_yr) + + # Pass #2: Develop current countries to fill remaining demand from the pass 1 + pass2_use <- pass1 %>% + # Reduce to undeveloped cells + filter(developed=="no") %>% + # Add remaining production available per period + left_join(pass1_leftover %>% select(period, prod_mt_yr_fill), by="period") %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_yr_fill, "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge the two passes + data2 <- bind_rows(pass1_use, pass2_use) + + } + + } + + + # Optimum development + if(dev_scenario_do=="Optimum"){ + + + # Calculate meat/person from all sectors in 2017 + all2017 <- hist %>% + select(-prod_mt) %>% + rename(npeople2017=npeople, meat_mt2017=meat_mt, meat_kg_person2017=meat_kg_person_2017) + + # Calculate meat/person from capture fisheries in 2100 + wc2100 <- wc %>% + # Reduce to climate scenario + filter(rcp==rcp_do & scenario=="Full Adaptation" & year %in% 2091:2100) %>% + # Calculate production, meat, meat/person + group_by(country, iso3) %>% + summarize(npeople2100=mean(npeople, na.rm=T), + meat_mt2100=mean(meat_mt, na.rm=T), + meat_kg_person2100=mean(meat_kg_person, na.rm=T)) + + # Merge + props <- all2017 %>% + left_join(wc2100) %>% + # Calculate missing meat + mutate(meat_mt_need=meat_kg_person2017/1000*npeople2100, + meat_mt_miss=meat_mt_need - meat_mt2100) %>% + # Remove countries gaining seafood/capita + filter(meat_mt_miss>0 * !is.na(meat_mt_miss)) %>% + # Remove countries without aquaculture potential + filter(iso3 %in% unique(data_orig$ter1_iso)) %>% + # Arrange in order of ease + arrange(meat_mt_miss) + + # Old way Calculate meat/person from capture fisheries in 2100 + # props <- wc %>% + # # Reduce to climate scenario + # filter(rcp==rcp_do & scenario=="Full Adaptation") %>% + # # Calculate amount of seafood needed to break even + # group_by(country, iso3) %>% + # summarize(meat_mt_need = meat_kg_person[year==2012]/1000*npeople[year==2100], + # meat_mt_avail = meat_mt[year==2100], + # meat_mt_miss = meat_mt_need - meat_mt_avail) %>% + # ungroup() %>% + # # Remove countries gaining seafood/capita + # filter(meat_mt_miss>0 * !is.na(meat_mt_miss)) %>% + # # Remove countries without aquaculture potential + # filter(iso3 %in% unique(data_orig$ter1_iso)) %>% + # # Arrange in order of ease + # arrange(meat_mt_miss) + + # Loop through countries in order of ease to fill + data2 <- purrr::map_df(1:nrow(props), function(i){ + + # Country do + iso3_do <- props$iso3[i] + + # Develop cells until satisfied + cdata <- data1 %>% + # Calculate meat + mutate(meat_mt_yr=0.87*prod_mt_yr) %>% + # Add meat needed + left_join(select(props, iso3, meat_mt_miss), by=c("ter1_iso_use"="iso3")) %>% + filter(!is.na(meat_mt_miss)) %>% + # Reduce to country of interest + filter(ter1_iso_use %in% iso3_do) %>% + # Arrange by profitability + group_by(period) %>% + arrange(period, desc(profits_usd_yr)) %>% + # Calculate cumulative production and see when it exceeds need + mutate(cum_meat_mt_yr=cumsum(meat_mt_yr), + over=cum_meat_mt_yr>meat_mt_miss, + ncells_over=cumsum(over)) %>% + # Reduce to only cells necessary + filter(ncells_over<=1) %>% + mutate(developed="yes") %>% + ungroup() %>% + # Add rank + mutate(dev_rank=i, + dev_order=1:n()) + + }) + + # Restrict by feed then restrict by demand + pass1 <- data2 %>% + # Arrange by period and development order (easy countries first) + arrange(period, dev_rank, dev_order) %>% + # Calculate cumulative feed and see if there is enough feed + group_by(period) %>% + mutate(cum_ff_mt_yr=cumsum(ff_demand_mt_yr)) %>% + mutate(feed_avail=ifelse(cum_ff_mt_yr<=ff_avail_mt_yr, "yes", "no")) %>% + # Calculate cumulative production and see if there is enough demand + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + demand_avail=ifelse(cum_prod_mt_yr<=prod_mt_cap, "yes", "no")) %>% + # Developed + mutate(developed=ifelse(demand_avail=="yes" & feed_avail=="yes", "yes", "no")) %>% + filter(developed=="yes") + + # Develop other countries to fill gap, if any + avail_key <- pass1 %>% + group_by(period) %>% + summarize(ff_mt_used=sum(ff_demand_mt_yr), + prod_mt_used=sum(prod_mt_yr)) + + # Pass 2 + pass2 <- data1 %>% + # Remove countries already focused on + filter(!ter1_iso_use %in% props$iso3) %>% + # Add amount of feed and production that is still available + left_join(avail_key) %>% + mutate(ff_mt_left=ff_avail_mt_yr - ff_mt_used, + prod_mt_left=prod_mt_cap - prod_mt_used) %>% + # Group by period, eez name + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + # Calculate cumulative feed and demand + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + cum_ff_mt_yr=cumsum(ff_demand_mt_yr), + ff_avail=ifelse(cum_ff_mt_yr <= ff_mt_left, "yes", "no"), + demand_avail=ifelse(cum_prod_mt_yr <= prod_mt_left, "yes", "no")) %>% + # Classify develop + mutate(developed=ifelse(ff_avail=="yes" & demand_avail == "yes", "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge pass 1 and pass 2 + data3 <- bind_rows(pass1, pass2) + data2 <- data3 + + } # end optimum + + # Perform final formatting + ############################################################### + + # Format data + data3 <- data2 %>% + # Add columns + mutate(mgmt_scenario=mgmt_scenario, + feed_scenario=feed_scenario, + dev_scenario=dev_scenario) %>% + # Arrange columns + # BIG NOTE: you get rid of a lot of columns here + select(rcp, mgmt_scenario, feed_scenario, dev_scenario, period, + sov1_iso, sov1_name, ter1_continent, ter1_iso_use, ter1_name_use, x, y, + species, feed_group, fifo_use, + prod_mt_yr:profits_usd_yr, ff_demand_mt_yr) %>% + rename(ter1_iso=ter1_iso_use, ter1_name=ter1_name_use) + + # Plot results + ############################################################### + + # Calculate statistics + stats <- data3 %>% + group_by(period, ter1_continent) %>% + summarize(ff_mmt_yr=sum(ff_demand_mt_yr)/1e6, + prod_mmt_yr=sum(prod_mt_yr)/1e6) + + # Plot production statistics + scen_title <- paste(rcp_do, mgmt_scenario_do, feed_scenario_do, dev_scenario_do, sep=", ") + g <- ggplot(stats, aes(x=period, y=prod_mmt_yr, fill=ter1_continent)) + + geom_bar(stat="identity") + + labs(x="", y="Production (millions mt)", title=scen_title) + + geom_hline(yintercept=prod_mt_cap/1e6) + + theme_bw() + print(g) + + # Export data + return(data3) + +} + + +# Develop demand-limited finfish mariculture +################################################################################ + + + +# Do all scenarios - new costs +##################################### + +# Costs scalars +cost_scalars <- c(2.1, 2.3, 2.5, 2.7) + + +# Build scenario key +rcps <- paste("RCP", c("2.6", "4.5", "6.0", "8.5")) +mgmt_scens <- c("BAU fisheries management", "Reformed fisheries management") +feed_scens <- c("BAU feed use", "Reformed feed use") +dev_scens <- c("Current", "Proportional", "Need-based") +# dev_scens <- c("Optimum") +scen_key <- expand.grid(rcp=rcps, + mgmt_scenario=mgmt_scens, + feed_scenario=feed_scens, + dev_scenario=dev_scens) %>% + arrange(rcp, mgmt_scenario, feed_scenario, dev_scenario) + +# Loop through cost scalars +for(i in 1:length(cost_scalars)){ + + # Suffix + suffix <- paste0("cost_search_", cost_scalars[i]) + + # Develop mariculture + output <- purrr::map_df(1:nrow(scen_key), function(x) { + + # Scenario + rcp_do <- scen_key$rcp[x] + mgmt_do <- scen_key$mgmt_scenario[x] + feed_do <- scen_key$feed_scenario[x] + dev_do <- scen_key$dev_scenario[x] + + # Expand mariculture + out_df <- expand_mariculture(rcp=rcp_do, mgmt_scenario = mgmt_do, feed_scenario = feed_do, dev_scenario=dev_do, suffix=suffix) + + }) + + # Export output + outfile <- paste0("finfish_output_", suffix, ".Rds") + saveRDS(output, file=file.path(outputdir, outfile)) + +} + diff --git a/code/cost_search/Step5_build_bivalve_scenario_results.R b/code/cost_search/Step5_build_bivalve_scenario_results.R new file mode 100644 index 0000000..1c88324 --- /dev/null +++ b/code/cost_search/Step5_build_bivalve_scenario_results.R @@ -0,0 +1,423 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(raster) +library(ggplot2) +library(tidyverse) +library(countrycode) +library(rnaturalearth) + +# Directories +datadir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/processed" +outputdir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_processed" +plotdir <- "figures" +tabledir <- "tables" + +# Read species data +load("data/species/data/aquaculture_species_key.Rdata") +spp_key <- data +rm(data_full, data) + +# Read feed info +maq_prod <- read.csv("data/feed_params/processed/FAO_2013_2017_maq_prod_averages_by_country.csv", as.is=T) + +# World +world <- rnaturalearth::ne_countries(scale="large", type = "countries", returnclass = "sf") + +# Read EEZ key +eezdir <- "data/eezs" +eez_key <- read.csv(file.path(eezdir, "eezs_v10_key.csv"), as.is=T) + +# Read population data +popdir <- "data/pop_growth/data" +pop <- readRDS(file.path(popdir, "WB_UN_1960_2100_human_population_by_country.Rds")) + +# Read fisheries impacts +wc <- readRDS("data/capture_projections/data/Free_etal_2020_national_projections_with_pop_data.Rds") + + +# Function to develop demand-limited bivalve mariculture +################################################################################ + +# For testing +rcp <- "RCP 2.6" +dev_scenario <- "Need-based" +suffix <- "cost_search_1.3" + +# Function to develop feed-constrained finfish mariculture +expand_mariculture <- function(rcp, dev_scenario, suffix=""){ + + # Steps + # 1. Determine feed availability + # 2. Determine feed demand + # 3. Develop based on distribution pattern + + # Setup + ############################################################### + + # Parameters + rcp_do <- rcp + dev_scenario_do <- dev_scenario + + # Demand limit + prod_mt_cap <- 103 * 0.17 / 0.17 * 1e6 + + # Read finfish forecasts + rcp_do_short <- gsub("\\.| ", "", rcp_do) + if(suffix==""){ + infile <- paste0( rcp_do_short, "_Bivalve_rational_use.Rds") + }else{ + infile <- paste0( rcp_do_short, "_Bivalve_rational_use_", suffix, ".Rds") + } + data_orig <- readRDS(file.path(outputdir, infile)) + + # Build forecasts period key + periods <- data_orig %>% pull(period) %>% unique() + period_key <- purrr::map_df(periods, function(x){ + yr1 <- substr(x,1,4) %>% as.numeric() + yr2 <- substr(x,6,9) %>% as.numeric() + yrs <- yr1:yr2 + yr_df <- tibble(period=x, + year=yrs) + }) + + + + # Develop mariculture based on development scenario + ############################################################### + + # Rational development + if(dev_scenario_do=="Rational"){ + + # Develop profitable cells first + data1 <- data_orig %>% + group_by(period) %>% + arrange(period, desc(profits_usd_yr)) %>% + # Calculate cumulative production + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_cap, "yes", "no")) %>% + filter(developed=="yes") + + } + + # Current development + if(dev_scenario_do=="Current"){ + + # Format current producers + # Based on FAO 2013-2017 averages + props <- maq_prod %>% + # Reduce to finfish aquaculture totals + filter(major_group=="Mollusca") %>% + # Format country information + mutate(country_use=countrycode(country_use, "country.name", "country.name"), + iso3_use=countrycode(country_use, "country.name", "iso3c")) %>% + # Remove TERRITORIES not forecast to have to viable mariculutre + filter(iso3_use %in% unique(data_orig$ter1_iso)) %>% + # Calculate proportion of feed dedicated to this TERRITORY + mutate(prod_prop=prod_mt/sum(prod_mt)) # note if a country isn't present in a period some feed will go unused + + # Pass #1. Develop profitable cells first, in countries with production + pass1 <- data_orig %>% + # Add percent of production allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3_use, prod_prop), by=c("ter1_iso"="iso3_use")) %>% + # Calculate production allocation + mutate(prod_prop=ifelse(is.na(prod_prop), 0, prod_prop), + prod_allocated_mt_yr=prod_mt_cap*prod_prop) %>% + # Remove places not allocated any production + filter(prod_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name) %>% + arrange(period, ter1_name, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative production and see if it exceeds ALLOCATED AMOUNT + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_allocated_mt_yr, "yes", "no")) %>% + ungroup() + + # Cells developed in pass one + pass1_use <- pass1 %>% + filter(developed=="yes") + + # Undeveloped production from pass 1 + pass1_leftover <- pass1_use %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) %>% + ungroup() %>% + mutate(prod_mt_yr_fill=prod_mt_cap-prod_mt_yr) + + # Pass #2: Develop current countries to fill remaining demand from the pass 1 + pass2_use <- pass1 %>% + # Reduce to undeveloped cells + filter(developed=="no") %>% + # Add remaining production available per period + left_join(pass1_leftover %>% select(period, prod_mt_yr_fill), by="period") %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_yr_fill, "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge the two passes + data1 <- bind_rows(pass1_use, pass2_use) + + } + + # Proportional development + if(dev_scenario_do=="Proportional"){ + + # 2050 population + pop55 <- pop %>% + filter(year==2100) + + # ISO key for production forecasts + props <- data_orig %>% + select(sov1_iso, sov1_name, ter1_iso, ter1_name) %>% + unique() %>% + arrange(ter1_iso) %>% + # Add 2055 population estimates + left_join(pop55 %>% select(iso3, pop_size_50perc), by=c("ter1_iso"="iso3")) %>% + # Only consider areas with population projection estimates + filter(!is.na(pop_size_50perc)) %>% + # Calculate proportional allocation + mutate(prod_prop=pop_size_50perc/sum(pop_size_50perc)) %>% + # Simplify + select(ter1_iso, prod_prop) %>% + # Rename ISO column so that I can use the same code as in the CURRENT scenario + rename(iso3_use=ter1_iso) + + # Pass #1. Develop profitable cells first, in countries with production + pass1 <- data_orig %>% + # Add percent of production allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3_use, prod_prop), by=c("ter1_iso"="iso3_use")) %>% + # Calculate production allocation + mutate(prod_prop=ifelse(is.na(prod_prop), 0, prod_prop), + prod_allocated_mt_yr=prod_mt_cap*prod_prop) %>% + # Remove places not allocated any production + filter(prod_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name) %>% + arrange(period, ter1_name, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative production and see if it exceeds ALLOCATED AMOUNT + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_allocated_mt_yr, "yes", "no")) %>% + ungroup() + + # Cells developed in pass one + pass1_use <- pass1 %>% + filter(developed=="yes") + + # Undeveloped production from pass 1 + pass1_leftover <- pass1_use %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) %>% + ungroup() %>% + mutate(prod_mt_yr_fill=prod_mt_cap-prod_mt_yr) + + # Pass #2: Develop current countries to fill remaining demand from the pass 1 + pass2_use <- pass1 %>% + # Reduce to undeveloped cells + filter(developed=="no") %>% + # Add remaining production available per period + left_join(pass1_leftover %>% select(period, prod_mt_yr_fill), by="period") %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_yr_fill, "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge the two passes + data1 <- bind_rows(pass1_use, pass2_use) + + } + + # Need-based development + if(dev_scenario_do=="Need-based"){ + + # Build proportions + props <- wc %>% + # Reduce to climate scenario + filter(rcp==rcp_do & scenario=="Full Adaptation") %>% + # Calculate amount of seafood needed to break even + group_by(country, iso3) %>% + summarize(meat_mt_need = meat_kg_person[year==2012]/1000*npeople[year==2100], + meat_mt_avail = meat_mt[year==2100], + meat_mt_miss = meat_mt_need - meat_mt_avail) %>% + ungroup() %>% + # Remove countries gaining seafood/capita + filter(meat_mt_miss>0 * !is.na(meat_mt_miss)) %>% + # Remove countries without aquaculture potential + filter(iso3 %in% unique(data_orig$ter1_iso)) %>% + # Calculate proportion of feed dedicated to this TERRITORY + mutate(prod_prop=meat_mt_miss/sum(meat_mt_miss), + prod_prop=1/n()) %>% + # Rename ISO column so that I can use the same code as in the CURRENT scenario + rename(iso3_use=iso3) + + # Pass #1. Develop profitable cells first, in countries with production + pass1 <- data_orig %>% + # Add percent of production allocated to TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + left_join(select(props, iso3_use, prod_prop), by=c("ter1_iso"="iso3_use")) %>% + # Calculate production allocation + mutate(prod_prop=ifelse(is.na(prod_prop), 0, prod_prop), + prod_allocated_mt_yr=prod_mt_cap*prod_prop) %>% + # Remove places not allocated any production + filter(prod_prop>0) %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) + group_by(period, ter1_name) %>% + arrange(period, ter1_name, rank_in_eez, desc(profits_usd_yr)) %>% + # Calculate cumulative production and see if it exceeds ALLOCATED AMOUNT + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_allocated_mt_yr, "yes", "no")) %>% + ungroup() + + # Cells developed in pass one + pass1_use <- pass1 %>% + filter(developed=="yes") + + # Undeveloped production from pass 1 + pass1_leftover <- pass1_use %>% + group_by(period) %>% + summarize(prod_mt_yr=sum(prod_mt_yr)) %>% + ungroup() %>% + mutate(prod_mt_yr_fill=prod_mt_cap-prod_mt_yr) + + # Pass #2: Develop current countries to fill remaining demand from the pass 1 + pass2_use <- pass1 %>% + # Reduce to undeveloped cells + filter(developed=="no") %>% + # Add remaining production available per period + left_join(pass1_leftover %>% select(period, prod_mt_yr_fill), by="period") %>% + # Rank by profitability within EEZ + group_by(period, eez_name) %>% + arrange(period, eez_name, desc(profits_usd_yr)) %>% + mutate(rank_in_eez=1:n()) %>% + ungroup() %>% + # Arrange within TERRITORY-LIKE UNIT (USA (Alaska/Hawaii/Continental), Guam, Puerto Rico all seperate) - basically distribute to current producers equally + group_by(period) %>% + arrange(period, rank_in_eez) %>% + mutate(cum_prod_mt_yr=cumsum(prod_mt_yr), + developed=ifelse(cum_prod_mt_yr <= prod_mt_yr_fill, "yes", "no")) %>% + ungroup() %>% + # Only developed cells + filter(developed=="yes") + + # Merge the two passes + data1 <- bind_rows(pass1_use, pass2_use) + + + + } + + # Perform final formatting + ############################################################### + + # Format data + data2 <- data1 %>% + # Add columns + mutate(rcp=rcp_do, + dev_scenario=dev_scenario_do) %>% + # Arrange columns + # BIG NOTE: you get rid of a lot of columns here + select(rcp, dev_scenario, period, + sov1_iso, sov1_name, ter1_continent, ter1_iso, ter1_name, x, y, + species, + prod_mt_yr:profits_usd_yr) + + # Plot results + ############################################################### + + # Calculate statistics + stats <- data2 %>% + group_by(period, ter1_continent) %>% + summarize(prod_mmt_yr=sum(prod_mt_yr)/1e6) + + # Plot production statistics + scen_title <- paste(rcp_do, dev_scenario_do, sep=", ") + g <- ggplot(stats, aes(x=period, y=prod_mmt_yr, fill=ter1_continent)) + + geom_bar(stat="identity") + + labs(x="", y="Production (millions mt)", title=scen_title) + + geom_hline(yintercept=prod_mt_cap/1e6, linetype="dashed") + + theme_bw() + print(g) + + # Export data + return(data2) + +} + + +# Develop demand-limited bivalve mariculture +################################################################################ + +# Costs scalars +cost_scalars <- c(2.1, 2.3, 2.5, 2.7) + +# Build scenario key +rcps <- paste("RCP", c("2.6", "4.5", "6.0", "8.5")) +dev_scens <- c("Current", "Proportional", "Need-based") +scen_key <- expand.grid(rcp=rcps, + dev_scenario=dev_scens) %>% + arrange(rcp, dev_scenario) + + +# Loop through cost scalars +for(i in 1:length(cost_scalars)){ + + # Suffix + suffix <- paste0("cost_search_", cost_scalars[i]) + + # Develop mariculture + output <- purrr::map_df(1:nrow(scen_key), function(x) { + + # Scenario + rcp_do <- scen_key$rcp[x] + dev_do <- scen_key$dev_scenario[x] + + # Expand mariculture + out_df <- expand_mariculture(rcp=rcp_do, dev_scenario=dev_do, suffix=suffix) + + }) + + # Export output + outfile <- paste0("bivalve_output_", suffix, ".Rds") + saveRDS(output, file=file.path(outputdir, outfile)) + +} + + + diff --git a/code/cost_search/Step6_merge_capture_mariculture_results.R b/code/cost_search/Step6_merge_capture_mariculture_results.R new file mode 100644 index 0000000..2ddbf88 --- /dev/null +++ b/code/cost_search/Step6_merge_capture_mariculture_results.R @@ -0,0 +1,145 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(raster) +library(ggplot2) +library(tidyverse) +library(countrycode) +library(rnaturalearth) + +# Directories +datadir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_processed" +plotdir <- "figures" +tabledir <- "tables" + +# Read capture projections +wc_orig <- readRDS("data/capture_projections/data/Free_etal_2020_global_projections_with_pop_data.Rds") + + +# Cost scalars +################################################################################ + +# Cost scalars +cost_scalars <- c(2.1, 2.3, 2.5, 2.7) + +# Loop through catch scalars +for(i in 1:length(cost_scalars)){ + + # Suffix + suffix <- paste0("cost_search_", cost_scalars[i]) + + # Read mariculture projections + faq_orig <- readRDS(file.path(datadir, paste0("finfish_output_", suffix, ".Rds"))) + baq_orig <- readRDS(file.path(datadir, paste0("bivalve_output_", suffix, ".Rds"))) + + # Global estimates + ################################################################################ + + # Period key + period_key <- tibble(year=c(2021:2030, 2051:2060, 2091:2100), + period=sort(rep(c("2021-2030", "2051-2060", "2091-2100"), 10))) + + # Calculate population growth + pop_global <- wc_orig %>% + # Select + select(year, npeople, npeople05, npeople20, npeople80, npeople95) %>% + # Add period + filter(year%in%c(2020, 2025, 2030, + 2050, 2055, 2060, + 2090, 2095, 2100)) %>% + unique() %>% + # Add period and summarize by period + mutate(period=ifelse(year<=2030, "2021-2030", + ifelse(year<=2060, "2051-2060", "2091-2100"))) %>% + group_by(period) %>% + summarize(npeople=mean(npeople), + npeople05=mean(npeople05), + npeople20=mean(npeople20), + npeople80=mean(npeople80), + npeople95=mean(npeople95)) + + # Calculate global statistics: capture + wc_global <- wc_orig %>% + # Add period + left_join(period_key) %>% + filter(!is.na(period)) %>% + # Summarize by period + group_by(rcp, scenario, period) %>% + summarize(prod_mt=mean(catch_mt), + meat_mt=mean(meat_mt)) %>% + ungroup %>% + # Simplify + mutate(sector="Capture fisheries", + scenario=recode(scenario, + "No Adaptation"="Business-as-usual", + "Full Adaptation"="Progressive reforms")) %>% + select(rcp, scenario, period, sector, prod_mt, meat_mt) + + # Calculate global statistics: finfish + faq_global <- faq_orig %>% + # Calculate global production + group_by(rcp, mgmt_scenario, feed_scenario, dev_scenario, period) %>% + summarise(prod_mt=sum(prod_mt_yr)) %>% + ungroup() %>% + # Calculate meat production + mutate(meat_mt=prod_mt*0.82) %>% + # Mark big-picture scenarios + mutate(scenario=ifelse(mgmt_scenario=="BAU fisheries management" & feed_scenario == "BAU feed use" & dev_scenario == "Current", "Business-as-usual", NA), + scenario=ifelse(mgmt_scenario=="Reformed fisheries management" & feed_scenario == "Reformed feed use" & dev_scenario == "Proportional", "Progressive reforms", scenario)) %>% + # Reduce to BAU or progressive reform + filter(!is.na(scenario)) %>% + # Simplify + mutate(sector="Finfish mariculture") %>% + select(rcp, scenario, period, sector, prod_mt, meat_mt) + + # Calculate global statistics: bivalve + baq_global_bau <- baq_orig %>% + # Calculate global production + group_by(rcp, dev_scenario, period) %>% + summarise(prod_mt=sum(prod_mt_yr)) %>% + ungroup() %>% + # Calculate meat production + mutate(meat_mt=prod_mt*0.17) %>% + # Reduce to current development patterns + filter(dev_scenario=="Current") %>% + # Simplify + mutate(sector="Bivalve mariculture", + scenario="Business-as-usual") %>% + select(rcp, scenario, period, sector, prod_mt, meat_mt) + baq_global_ref <- baq_global_bau %>% + mutate(scenario="Progressive reforms") + baq_global <- bind_rows(baq_global_bau, baq_global_ref) + + # Merge sectors + gdata <- bind_rows(wc_global, faq_global, baq_global) %>% + # Add population growth + left_join(pop_global, by="period") %>% + # Calculate meat per capita + mutate(meat_kg_person=meat_mt*1000/npeople, + meat_kg_person05=meat_mt*1000/npeople05, + meat_kg_person20=meat_mt*1000/npeople20, + meat_kg_person80=meat_mt*1000/npeople80, + meat_kg_person95=meat_mt*1000/npeople95) + + # Export + outfile <- paste0("global_capture_mariculture_output_merged_", suffix, ".Rds") + saveRDS(gdata, file=file.path(datadir, outfile)) + + + # Plot data + ################################################################################ + + g <- ggplot(gdata, aes(x=period, y=meat_kg_person, fill=sector)) + + facet_grid(rcp~scenario) + + geom_bar(stat="identity") + + labs(x="Period", y="Edible meat per capita (kg/person)") + g + + +} + diff --git a/code/cost_search/figure_cost_search_analysis.R b/code/cost_search/figure_cost_search_analysis.R new file mode 100644 index 0000000..74f6e37 --- /dev/null +++ b/code/cost_search/figure_cost_search_analysis.R @@ -0,0 +1,72 @@ + +# Clear workspace +rm(list = ls()) + +# Setup +################################################################################ + +# Packages +library(tidyverse) + +# Directories +datadir <- "/Volumes/GoogleDrive/Shared drives/emlab/projects/current-projects/blue-paper-2/data/output/cost_search_processed" +plotdir <- "figures" +tabledir <- "tables" + + +# Build data +################################################################################ + +# Cost scalars +cost_scalars <- c(1.3, 1.5, 1.7, 2, 2.1) + +# Loop through cost scalars and build data +x <- 1.3 +data_orig <- purrr::map_df(cost_scalars, function(x){ + + # Read data + infile <- paste0("global_capture_mariculture_output_merged_cost_search_", x, ".Rds") + fdata_orig <- readRDS(file.path(datadir, infile)) + + # Summarize data + fdata <- fdata_orig %>% + group_by(rcp, scenario, period) %>% + summarize(meat_kg_person=sum(meat_kg_person)) %>% + ungroup() %>% + mutate(cost_scalar=x) %>% + select(cost_scalar, everything()) + +}) + +# Data +data <- data_orig %>% + filter(scenario=="Progressive reforms" & period=="2091-2100") + + +# Plot data +################################################################################ + +g <- ggplot(data, aes(x=cost_scalar, y=meat_kg_person, color=rcp)) + + geom_line() + + geom_point() + + # Reference line + geom_hline(yintercept=8.918661) + + # Labels + labs(x="Cost scalar", y="Seafood production per capita") + + # Axis + scale_y_continuous(lim=c(0,NA)) + + # Theme + theme_bw() +g + + + + + + + + + + + +