1 Read in the raw data

This includes the FAO capture production data and a list of the “target” species.

## FAO capture production data - all columns being parsed as characters and producing error in one column, but not sure which? (read.csv might help avoid this error?)
fis_fao_raw <-  read_csv(file.path(dir_M, 'git-annex/globalprep/_raw_data/FAO_capture/d2019/Global_capture_production_Quantity_1950-2017.csv'))

# species list - used same raw files from v2018
sp2grp <-  read_csv(here('globalprep/prs_targetedharvest/v2019/raw/species2group.csv')) %>%
  dplyr::filter(incl_excl == 'include') %>%
  dplyr::select(target, species); head(sp2grp)

2 Clean the FAO data

# Rename columns and remove unit column
fao_clean <- fis_fao_raw %>% 
  dplyr::rename(country = "Country (Country)",
                species = "Species (ASFIS species)",
                area = "Fishing area (FAO major fishing area)") %>%
  select(-"Unit (Unit)")

# Gather by year and value to expand and make each line a single observation for country, species and year (tidy data!) 
fao_clean <- fao_clean %>%
  tidyr::gather("year", "value", -(1:3)) %>%
  dplyr::mutate(year = gsub("X", "", year)) %>%
    fao_clean_data() 

fao_clean <- fao_clean %>%
  dplyr::mutate(species = as.character(species)) %>%
  dplyr::mutate(species = ifelse(stringr::str_detect(species, "Henslow.*s swimming crab"), "Henslow's swimming crab", species))

3 Identify the target species

This analysis only includes target species. The warning messages need to be checked and, if necessary, changes should be made to the raw/species2group.csv

# check for discrepancies in species list
spgroups <-  sort(as.character(unique(fao_clean$species))) # species groups in FAO data 
groups <-  c('turtle', 'whale', 'dolphin', 'porpoise') # seals and sea lions removed from vector (pinnipeds no longer included) 

# Going through FAO data species and seeing if they're in our master list of species
## Looking to see if we need to add species that have changed name
### v2019: added Andrew's beaked whale to species2group.csv based on error message. Other species in error message are excluded.
for (group in groups) {# group='dolphin'
possibles <- spgroups[grep(group, spgroups)]
d_missing_l <-  setdiff(possibles, sp2grp$species)
  if (length(d_missing_l)>0){
    cat(sprintf("\nMISSING in the lookup the following species in target='%s'.\n    %s\n", 
                group, paste(d_missing_l, collapse='\n    ')))
  }
}


# check for species in lookup not found in data
l_missing_d <-  setdiff(sp2grp$species, spgroups)
if (length(l_missing_d)>0){
  cat(sprintf('\nMISSING: These species in the lookup are not found in the FAO data \n'))
  print(l_missing_d)
}


## filter data to include only target species ----
target_spp <-  fao_clean %>%
  dplyr::filter(species %in% sp2grp$species) # this goes from 2144 spp in FAO list to just 69

unique(target_spp$area) # confirm these are all marine regions

4 Summarize data

# widen spread to expand years
wide = target_spp %>%
  tidyr::spread(year, value) %>%
  dplyr::left_join(sp2grp, by='species'); head(wide)


# gather long by target
long = wide %>%
  dplyr::select(-area) %>%
  tidyr::gather(year, value, -country, -species, -target, na.rm=T) %>%
  dplyr::mutate(year = as.integer(as.character(year))) %>%
  dplyr::arrange(country, target, year); head(long)


# explore Japan[210] as an example
japan <- long %>% 
  dplyr::group_by(country, target, year) %>%
  dplyr::summarize(value = sum(value)) %>% 
  dplyr::filter(country == 'Japan', target == 'cetacean', year >= 2000) 

# summarize totals per region per year - number of individual animals from each spp group? 
sum = long %>%
  dplyr::group_by(country, year) %>%
  dplyr::summarize(value = sum(value, na.rm=TRUE)) %>%
  dplyr::filter(value != 0) %>%
  dplyr::ungroup(); head(sum) 

5 Assign country names to OHI regions

sum <- sum %>%
  dplyr::mutate(country = as.character(country)) %>%
  dplyr::mutate(country = ifelse(stringr::str_detect(country, "C.*te d'Ivoire"), "Ivory Coast", country))


### Function to convert to OHI region ID
m_sum_rgn <- name_2_rgn(df_in = sum, 
                       fld_name='country', 
                       flds_unique=c('year'))

# Filter out duplicates based on error message from previous step
dplyr::filter(m_sum_rgn, country %in% c("Guadeloupe", "Martinique"))

# They will be summed:
m_sum_rgn <- m_sum_rgn %>%
  dplyr::group_by(rgn_id, rgn_name, year) %>%
  dplyr::summarize(value = sum(value)) %>%
  dplyr::ungroup()

6 Scale the data and save files

Data is rescaled by dividing by the 95th quantile of values across all regions from 2011 to 2017 (most recent year of FAO data).

target_harvest <- m_sum_rgn %>%
  dplyr::mutate(quant_95 = quantile(value[year %in% 2011:2017], 0.95, na.rm = TRUE)) %>%
  dplyr::mutate(score = value / quant_95) %>% 
  dplyr::mutate(score = ifelse(score>1, 1, score)) %>%
  dplyr::select(rgn_id, year, pressure_score = score) %>%
  dplyr::arrange(rgn_id, year); head(target_harvest); summary(target_harvest)
  
# v2019 quant_95 = 3407 


# any regions that did not have a catch should have score = 0 
rgns <-  rgn_master %>%
  dplyr::filter(rgn_typ == "eez") %>%
  dplyr::select(rgn_id = rgn_id_2013) %>%
  dplyr::filter(rgn_id < 255) %>%
  base::unique() %>%
  dplyr::arrange(rgn_id)

# This is just a list of rgn IDS - do we want to update it to a rgn list more recent than 2013?
  
rgns <- expand.grid(rgn_id = rgns$rgn_id, year = min(target_harvest$year):max(target_harvest$year))
  
target_harvest <-  rgns %>%
  dplyr::left_join(target_harvest) %>%
  dplyr::mutate(pressure_score = ifelse(is.na(pressure_score), 0, pressure_score)) %>%
  dplyr::arrange(rgn_id); head(target_harvest); summary(target_harvest)


write_csv(target_harvest, 
          file.path(here('globalprep/prs_targetedharvest/v2019/output/fao_targeted.csv')))


target_harvest_gf <- target_harvest %>%
  dplyr::mutate(gapfill = 0) %>%
  dplyr::select(rgn_id, year, gapfill)
# all zeroes for gapfill column; nothing being gapfilled but need to have a record 

write_csv(target_harvest_gf,
          file.path(here('globalprep/prs_targetedharvest/v2019/output/fao_targeted_gf.csv')))

7 Data check

The data from last year and this year should be the same unless there were changes to underlying FAO data or the master species list.

In this case, all of the regions looked very similar.

new <- read_csv(here("globalprep/prs_targetedharvest/v2019/output/fao_targeted.csv")) %>% 
  filter(year==2016)
# pull just 2015 data from target_harvest df - should we change this to a more recent year for v2019?

old <- read_csv(here("globalprep/prs_targetedharvest/v2018/output/fao_targeted.csv")) %>%
  dplyr::filter(year == 2016) %>%
  dplyr::select(rgn_id, year, pressure_score_old=pressure_score) %>%
  dplyr::left_join(new, by=c("rgn_id", "year"))


compare_plot <- ggplot(data = old, aes(x=pressure_score_old, y= pressure_score, label=rgn_id))+
  geom_point()+
  geom_abline(color="red")

plot(compare_plot)
ggplotly(compare_plot)