EWG-21-18 ToR3: update the economic scripts for the static analysis

Study 1: space-time analysis

rm(list=ls())
# path <- file.path("C:/Users/Richard.Curtin/OneDrive - bord iascaigh mhara/Documents/Dropbox Files/2014/BIM ECON/EU_Meetings/STECF EWG CelticSea 2021_2118/STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-1-CatchComposition")
# path_data <- file.path("C:/Users/Richard.Curtin/OneDrive - bord iascaigh mhara/Documents/Dropbox Files/2014/BIM ECON/EU_Meetings/STECF EWG CelticSea 2021_2118/STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-2-ThresholdsAnalysis")
user_path <- "C:/STECF_Celtic_Sea_EWG-21-18/Software/"
user_path_data <- "C:/STECF_Celtic_Sea_EWG-21-18/Software/"
path <- paste0(user_path, file.path("STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-1-CatchComposition"))
path_data <- paste0(user_path_data, file.path("STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-2-ThresholdsAnalysis"))

setwd(path)

load(file.path(path_data,"NWW_data_set_w_thresholds.RData"))

trips_in_cspz  <- unique(data[data$in_CSPZ==1,]$trip_id)
length(trips_in_cspz)
## [1] 35673
trips_out_cspz <- unique(data[data$in_CSPZ==0,]$trip_id)
length(trips_out_cspz)  # outside trips AND CSPZ overlapping trips
## [1] 297483
data$trip_w_overlap <- 0
data[(data$trip_id %in% trips_in_cspz) & (data$trip_id %in% trips_out_cspz), "trip_w_overlap"] <- 1
length(unique(data[data$trip_w_overlap==1, "trip_id"]))  # nb overlapping trips
## [1] 6749
trips_out_cspz <- trips_out_cspz[!trips_out_cspz %in% trips_in_cspz]
length(trips_out_cspz)  # entirely outside trips
## [1] 290734

Study 2: GLM statistical analysis

rm(list=ls())
# path <- file.path('C:/Users/Richard.Curtin/OneDrive - bord iascaigh mhara/Documents/Dropbox Files/2014/BIM ECON/EU_Meetings/STECF EWG CelticSea 2021_2118/STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-1-GLM-on-CPUEs')
# path_data <- file.path('C:/Users/Richard.Curtin/OneDrive - bord iascaigh mhara/Documents/Dropbox Files/2014/BIM ECON/EU_Meetings/STECF EWG CelticSea 2021_2118/STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-2-ThresholdsAnalysis')
user_path <- "C:/STECF_Celtic_Sea_EWG-21-18/Software/"
user_path_data <- "C:/STECF_Celtic_Sea_EWG-21-18/Software/"
path <- paste0(user_path, file.path("STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-1-GLM-on-CPUEs"))
path_data <- paste0(user_path_data, file.path("STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-2-ThresholdsAnalysis"))
setwd(path)

load(file.path(path_data,"NWW_data_set.RData"))

# Remind the list of rectangles in CSPZ
dd <- table(data$in_CSPZ, data$rect_cod)
rect_in <- colnames(dd[,dd[2,]!=0])
rect_out <- colnames(dd[,dd[1,]!=0])
rect_in # rectangle inside
##  [1] "28D9" "28E0" "28E1" "28E2" "29D9" "29E0" "29E1" "29E2" "29E3" "29E4"
## [11] "30D9" "30E0" "30E1" "30E2" "30E3" "30E4" "30E5" "31D9" "31E0" "31E1"
## [21] "31E2" "31E3" "31E4" "31E5" "32D9" "32E0" "32E1" "32E2" "32E3" "32E4"
## [31] "32E5" "33D9" "33E0"

Study 3: Threshold analysis

rm(list=ls())
# path <- file.path('C:/Users/Richard.Curtin/OneDrive - bord iascaigh mhara/Documents/Dropbox Files/2014/BIM ECON/EU_Meetings/STECF EWG CelticSea 2021_2118/STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-2-ThresholdsAnalysis')
user_path <- "C:/STECF_Celtic_Sea_EWG-21-18/Software/"
path <- paste0(user_path, file.path("STECF-EWG-21-18/ToR1/NWW-MS-dataset/ToR1-2-ThresholdsAnalysis"))
setwd(path)


# read and format
# Import data:
if(FALSE){
  data<-read.csv2('data_set_2.csv', sep=",")
  # format:
  an <- function(x) as.numeric(as.character(x))
  data$catch_value_euro      <- an(data$catch_value_euro)
  data$catch_volume_kg       <- an(data$catch_volume_kg)
  data$fishing_hours         <- an(data$fishing_hours)
  data$nb_people_employed    <- an(data$nb_people_employed)
  data$in_CSPZ               <- an(data$in_CSPZ)
  save(data, file='NWW_data_set.RData') # decrease the memory load by using an RData
}
load(file='NWW_data_set.RData')

#file.exists('NWW_data_set.RData')
#file.exists('data_set_2.csv')

#data<-read.csv2('data_set_2.csv', sep=",")
#head(data)

print(table(data$MS_cod, data$year))
##      
##          2017    2018    2019    2020
##   BEL   91887  100185   97391       0
##   ESP   33918   37560   36022   34607
##   FRA 1047397 1012412  961230       0
##   IRL  444777  412629  405714  350926
print(table(data$MS_cod=="FRA", data$gear_cod, data$year))
## , ,  = 2017
## 
##        
##                    BB    DRB    FDV    FOO    FPO    GEN    GES    GNC    GND
##   FALSE     19     93   5804      0      0   9989      0      0      0      8
##   TRUE      52      0  29128      2     10  60994      0     97     28      0
##        
##            GNE    GNS   GNS4    GTN    GTR    HMS     LH    LHM   LHM1    LHP
##   FALSE      0  15809      0      0    731      0      0    147      0     76
##   TRUE       0  40273      0    273  48864   1880     21   1217      0   5519
##        
##           LHP1   LHP2     LL    LLD    LLF    LLS   LLS1   LLS2    LNB    LTL
##   FALSE      0      0      0    542      0  11024      0      0      0     32
##   TRUE       0      0   1702     11     13  13168      0      0     23   2709
##        
##            LVS    MIS    OTB    OTM    OTT     PS    PTB    PTM    SDN    SDV
##   FALSE      0     40 151909   1202 187536     39     14   1915    543      0
##   TRUE       1      0 581515  13465 232611   2623      0   3098   5245   1737
##        
##            SPR    SSC     TB    TBB
##   FALSE      0  30151      7 152952
##   TRUE       0      0      0   1118
## 
## , ,  = 2018
## 
##        
##                    BB    DRB    FDV    FOO    FPO    GEN    GES    GNC    GND
##   FALSE     23      0   4617      0      0  10659      0      0      0     10
##   TRUE      75      0  28236      0      3  60353      0     98      0     13
##        
##            GNE    GNS   GNS4    GTN    GTR    HMS     LH    LHM   LHM1    LHP
##   FALSE      0  13046     91      0    896      0      0    103     30    230
##   TRUE       0  39294      0   1169  51908   2095     86   1103      0   5999
##        
##           LHP1   LHP2     LL    LLD    LLF    LLS   LLS1   LLS2    LNB    LTL
##   FALSE     19     53    195    388      0     56  12386    923      0    304
##   TRUE       0      0   1222      6      6   9398      0      0     49   2276
##        
##            LVS    MIS    OTB    OTM    OTT     PS    PTB    PTM    SDN    SDV
##   FALSE      0      0 128933    589 178742      0      0   1873      0      0
##   TRUE       0      0 566921  18269 214467   2245     67   3796    782   2080
##        
##            SPR    SSC     TB    TBB
##   FALSE      0  29594      0 166614
##   TRUE       0      0      0    396
## 
## , ,  = 2019
## 
##        
##                    BB    DRB    FDV    FOO    FPO    GEN    GES    GNC    GND
##   FALSE      2      0   4398      0      0  13760      0      0      0     76
##   TRUE      17      0  28734      1      4  58822     44    119    118      7
##        
##            GNE    GNS   GNS4    GTN    GTR    HMS     LH    LHM   LHM1    LHP
##   FALSE      0  12781    515     23    583      0      0     16     20    104
##   TRUE      73  39557      0   1278  54800   1853     16    783      0   6484
##        
##           LHP1   LHP2     LL    LLD    LLF    LLS   LLS1   LLS2    LNB    LTL
##   FALSE     10     21      0    159      0     35  12189    595      0    221
##   TRUE       0      0    605      4      2  10015      0      0     46   2151
##        
##            LVS    MIS    OTB    OTM    OTT     PS    PTB    PTM    SDN    SDV
##   FALSE      0      0 138246   1033 152436      0     25   1703   1049      0
##   TRUE       0      0 545160    845 200391   2167      0   3816    466   1781
##        
##            SPR    SSC     TB    TBB
##   FALSE      3  38949      0 160175
##   TRUE       0      0      0   1071
## 
## , ,  = 2020
## 
##        
##                    BB    DRB    FDV    FOO    FPO    GEN    GES    GNC    GND
##   FALSE      7      0   3888      0      0  12292      0      0      0    139
##   TRUE       0      0      0      0      0      0      0      0      0      0
##        
##            GNE    GNS   GNS4    GTN    GTR    HMS     LH    LHM   LHM1    LHP
##   FALSE      0  13828    244     77    318      0      0     76     27    301
##   TRUE       0      0      0      0      0      0      0      0      0      0
##        
##           LHP1   LHP2     LL    LLD    LLF    LLS   LLS1   LLS2    LNB    LTL
##   FALSE      8     56      0     74      0     10   9948    718      0    329
##   TRUE       0      0      0      0      0      0      0      0      0      0
##        
##            LVS    MIS    OTB    OTM    OTT     PS    PTB    PTM    SDN    SDV
##   FALSE      0      0 114687    990 128904      0    373   2589    332      0
##   TRUE       0      0      0      0      0      0      0      0      0      0
##        
##            SPR    SSC     TB    TBB
##   FALSE      0  32798      0  62520
##   TRUE       0      0      0      0
print(table(data$gear_cod, data$year, data$MS_cod=="FRA"))
## , ,  = FALSE
## 
##       
##          2017   2018   2019   2020
##            19     23      2      7
##   BB       93      0      0      0
##   DRB    5804   4617   4398   3888
##   FDV       0      0      0      0
##   FOO       0      0      0      0
##   FPO    9989  10659  13760  12292
##   GEN       0      0      0      0
##   GES       0      0      0      0
##   GNC       0      0      0      0
##   GND       8     10     76    139
##   GNE       0      0      0      0
##   GNS   15809  13046  12781  13828
##   GNS4      0     91    515    244
##   GTN       0      0     23     77
##   GTR     731    896    583    318
##   HMS       0      0      0      0
##   LH        0      0      0      0
##   LHM     147    103     16     76
##   LHM1      0     30     20     27
##   LHP      76    230    104    301
##   LHP1      0     19     10      8
##   LHP2      0     53     21     56
##   LL        0    195      0      0
##   LLD     542    388    159     74
##   LLF       0      0      0      0
##   LLS   11024     56     35     10
##   LLS1      0  12386  12189   9948
##   LLS2      0    923    595    718
##   LNB       0      0      0      0
##   LTL      32    304    221    329
##   LVS       0      0      0      0
##   MIS      40      0      0      0
##   OTB  151909 128933 138246 114687
##   OTM    1202    589   1033    990
##   OTT  187536 178742 152436 128904
##   PS       39      0      0      0
##   PTB      14      0     25    373
##   PTM    1915   1873   1703   2589
##   SDN     543      0   1049    332
##   SDV       0      0      0      0
##   SPR       0      0      3      0
##   SSC   30151  29594  38949  32798
##   TB        7      0      0      0
##   TBB  152952 166614 160175  62520
## 
## , ,  = TRUE
## 
##       
##          2017   2018   2019   2020
##            52     75     17      0
##   BB        0      0      0      0
##   DRB   29128  28236  28734      0
##   FDV       2      0      1      0
##   FOO      10      3      4      0
##   FPO   60994  60353  58822      0
##   GEN       0      0     44      0
##   GES      97     98    119      0
##   GNC      28      0    118      0
##   GND       0     13      7      0
##   GNE       0      0     73      0
##   GNS   40273  39294  39557      0
##   GNS4      0      0      0      0
##   GTN     273   1169   1278      0
##   GTR   48864  51908  54800      0
##   HMS    1880   2095   1853      0
##   LH       21     86     16      0
##   LHM    1217   1103    783      0
##   LHM1      0      0      0      0
##   LHP    5519   5999   6484      0
##   LHP1      0      0      0      0
##   LHP2      0      0      0      0
##   LL     1702   1222    605      0
##   LLD      11      6      4      0
##   LLF      13      6      2      0
##   LLS   13168   9398  10015      0
##   LLS1      0      0      0      0
##   LLS2      0      0      0      0
##   LNB      23     49     46      0
##   LTL    2709   2276   2151      0
##   LVS       1      0      0      0
##   MIS       0      0      0      0
##   OTB  581515 566921 545160      0
##   OTM   13465  18269    845      0
##   OTT  232611 214467 200391      0
##   PS     2623   2245   2167      0
##   PTB       0     67      0      0
##   PTM    3098   3796   3816      0
##   SDN    5245    782    466      0
##   SDV    1737   2080   1781      0
##   SPR       0      0      0      0
##   SSC       0      0      0      0
##   TB        0      0      0      0
##   TBB    1118    396   1071      0
sum(data$year==2019, data$MS_cod=="FRA")
## [1] 4521396
#in value in euro
round(rev(sort(tapply(data[which(data$inside_CSPZ=="Y"),"catch_value_euro"], data[which(data$inside_CSPZ=="Y"),"fao_cod"], sum)))) [1:15]
##      MNZ      WHG      CRW      HAD      COD      SOL      HER      MAC 
## 24126046 23420686 22181834 19899477 11319187  5694780  4080758  3641948 
##      SPR      RAJ      LBE      RJM      RJH      POK      PAL 
##  2899570  2331647  2307855  2297823  1677185  1421604  1229322
# in landed kg
round(rev(sort(tapply(data[which(data$inside_CSPZ=="Y"),"catch_volume_kg"], data[which(data$inside_CSPZ=="Y"),"fao_cod"], sum))))  [1:15]
##      HKE      WHG      HER      SPR      LEZ      NEP      ANF      HAD 
## 32770754 15409865 12827638 12759427 12707093 12235456 10811429 10629663 
##      MNZ      SCE      MAC      CRE      COD      WIT      ANE 
##  5445598  3896142  3849914  3783417  3592283  2898231  2774506
dd <- table(data$in_CSPZ, data$rect_cod)
rect_in <- colnames(dd[,dd[2,]!=0])
rect_out <- colnames(dd[,dd[1,]!=0])
head(data[data$rect_cod==" ",])   
##               X MS_cod vessel_id   trip_id year month nb_people_employed
## 3344431 3344431    ESP    ESP_v1 ESP_t1237 2018     8                 NA
## 3344432 3344432    ESP    ESP_v1 ESP_t1237 2018     8                 NA
## 3344433 3344433    ESP    ESP_v1 ESP_t1237 2018     8                 NA
## 3344461 3344461    ESP    ESP_v1 ESP_t1324 2018     9                 NA
## 3344462 3344462    ESP    ESP_v1 ESP_t1324 2018     9                 NA
## 3344474 3344474    ESP    ESP_v1 ESP_t1422 2018     9                 NA
##         gear_cod mesh_size selectivity_device metier_dcf_6_cod div_ciem_cod
## 3344431     LLS1        NA                 NA                      27.7.j.2
## 3344432     LLS1        NA                 NA                      27.7.j.2
## 3344433     LLS1        NA                 NA                      27.7.j.2
## 3344461     LLS1        NA                 NA                        27.7.h
## 3344462     LLS1        NA                 NA                        27.7.h
## 3344474     LLS1        NA                 NA                        27.7.h
##         rect_cod square square_coordinates inside_CSPZ fao_cod catch_volume_kg
## 3344431            <NA>                 NA           N     BLI           48.56
## 3344432            <NA>                 NA           N     BXD           47.00
## 3344433            <NA>                 NA           N     WRF            4.20
## 3344461            <NA>                 NA           N     ALF           11.10
## 3344462            <NA>                 NA           N     WRF            5.00
## 3344474            <NA>                 NA           N     BYS           40.50
##         catch_value_euro fishing_hours in_CSPZ
## 3344431            98.47            NA       0
## 3344432           752.05            NA       0
## 3344433            98.28            NA       0
## 3344461               NA            NA       0
## 3344462           117.00            NA       0
## 3344474           298.25            NA       0
tapply(data[data$rect_cod==" ",]$catch_volume_kg, data[data$rect_cod==" ",]$MS_cod, sum) # 1294 tons not spatially allocated from ESP, all gear types, but outside the CSPZ
##     BEL     ESP     FRA     IRL 
##      NA 1294650      NA      NA
tapply(data$catch_volume_kg, data$MS_cod, sum) # check per MS
##       BEL       ESP       FRA       IRL 
##  12351546 102493807 520970908 267131095

Study 4. Exploration

tables_plots_path <- "//galwayfs03/FishData/STECF/Meetings/STECF_technical_measures_2021/static_Richard/tables_plots_ToR3_StaticEconUpdate"
                    
sum(data[which(data$fao_cod=="COD"),]$catch_volume_kg)/1000
## [1] 4895.596
sum(data[data$fao_cod=="COD" & data$in_CSPZ==1,]$catch_volume_kg)/1000
## [1] 3713.513
length(unique(data$vessel_id))
## [1] 2047
length(unique(data$trip_id))
## [1] 326407
an <- function(x) as.numeric(as.character(x))
table(unlist(tapply(an(data$nb_people_employed), data$vessel_id, function (x) unique(is.na (x))))) # nb of vessels with informed nb crew (i.e. is.na at FALSE)
## 
## FALSE  TRUE 
##  1266   789
sum(tapply(an(data$nb_people_employed), data$vessel_id, function (x) max (x)), na.rm=TRUE)
## [1] 4463.53
trips_in_cspz  <- unique(data[data$in_CSPZ==1,]$trip_id)
length(trips_in_cspz)
## [1] 35673
trips_out_cspz <- unique(data[data$in_CSPZ==0,]$trip_id)
length(trips_out_cspz)  # outside trips AND CSPZ overlapping trips
## [1] 297483
data$trip_w_overlap <- 0
data[(data$trip_id %in% trips_in_cspz) & (data$trip_id %in% trips_out_cspz), "trip_w_overlap"] <- 1 
length(unique(data[data$trip_w_overlap==1, "trip_id"]))  # nb overlapping trips
## [1] 6749
trips_out_cspz <- trips_out_cspz[!trips_out_cspz %in% trips_in_cspz]  
length(trips_out_cspz)  # entirely outside trips
## [1] 290734
sum(data[data$fao_cod=="COD" & data$trip_w_overlap==1,]$catch_volume_kg)/1000
## [1] 1361.449
### missing code
dd <- tapply(data[data$fao_cod=="COD",]$catch_volume_kg, list(data[data$fao_cod=="COD",]$metier_dcf_6_cod), sum)
dd <- dd[order(dd, decreasing=TRUE)]
ddd <- cumsum(dd)/sum(dd, na.rm=TRUE)
nm <- names(ddd[ddd<0.99])
all_mets_with_cod <- nm[!is.na(nm)]
all_mets_with_cod
##  [1] "OTB_DEF_100_119_0"   "OTB_DEF_100-119_0_0" "OTT_DEF_100_119_0"  
##  [4] "SSC_DEF_100-119_0_0" "TBB_DEF_70-99_0_0"   "OTB_CRU_70-99_0_0"  
##  [7] "OTB_DEF_70_99_0"     "OTB_DEF_70-99_0_0"   "OTB_CRU_100-119_0_0"
## [10] "OTB_DEF_>=120_0_0"   "GNS_DEF_120-219_0_0" "OTT_CRU_100_119_0"  
## [13] "TBB_DEF_70-99"       "SDN_DEF_100_119_0"   "OTB_CEP_70_99_0"    
## [16] "OTM_DEF_100_119_0"   "OTB_MCD_70-99"       "GTR_DEF_>=220_0"
dd <- tapply(data[data$fao_cod=="COD",]$catch_volume_kg, list(data[data$fao_cod=="COD",]$metier_dcf_6_cod, data[data$fao_cod=="COD",]$year ), sum)
dd <- doBy::orderBy(~ - 2017, dd)
dd <- dd[all_mets_with_cod,]
all_OTB_OTT_PTB_mets_with_cod  <- all_mets_with_cod[grepl("OTB", all_mets_with_cod) | grepl("PTB", all_mets_with_cod) | grepl("OTT", all_mets_with_cod)]# | grepl("SSC", all_mets_with_cod) | grepl("SDN", all_mets_with_cod)]
# we remove Seines. Regulation is only for bottom trawlers

trips_in_cspz  <- unique(data[data$in_CSPZ==1 & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod,]$trip_id)

length(trips_in_cspz)
## [1] 16925

Economic loss if choke in 2019 for France

fra_cod <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="FRA" & fao_cod=="COD" & year==2019), sum)
fra_cod_sub <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="FRA" & fao_cod=="COD" & year == 2019 & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod), sum)
names(fra_cod) <- c("year", "catch_all_metiers")
names(fra_cod_sub) <- c("year", "all_OTB_OTT_PTB")
fra_cod_dif <- merge(fra_cod, fra_cod_sub)
fra_cod_dif$per <- fra_cod_dif$all_OTB_OTT_PTB/fra_cod_dif$catch_all_metiers
knitr::kable(fra_cod_dif)
year catch_all_metiers all_OTB_OTT_PTB per
2019 373250.4 355940.3 0.9536234
lgb1 <- data [data$year=="2019" & data$MS_cod=="FRA" & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod,]
lgb1 <- orderBy(~month, data=lgb1)
lgb1$is_cod     <- 0
lgb1$cumsum_cod <- 0
lgb1[lgb1$fao_cod=="COD", "is_cod"] <- 1
lgb1 [lgb1$is_cod==1, "cumsum_cod"] <- lgb1 [lgb1$is_cod==1, "catch_volume_kg"]
lgb1 [, "cumsum_cod"] <- cumsum(lgb1$cumsum_cod)
lgb1 [, "choked"] <- 0
lgb1 [lgb1 [, "cumsum_cod"] > 294000*fra_cod_dif$per, "choked"] <- 1
# lgb1 [lgb1 [, "cumsum_cod"] > 461000, "choked"] <- 1 # for IRL
weight_choked_not_choked <- tapply(lgb1$catch_volume_kg, list(lgb1$fao_cod, lgb1$choked), sum)
weight_choked_not_choked <- weight_choked_not_choked/1e3 # convert in tons
weight_choked_not_choked <- orderBy(~ -0, weight_choked_not_choked)
weight_choked_not_choked <- cbind.data.frame(weight_choked_not_choked,
                                             percentdiff= - weight_choked_not_choked[,"1"] /(weight_choked_not_choked[,"1"]+weight_choked_not_choked[,"0"])*100
)
euros_choked_not_choked  <- tapply(lgb1$catch_value_euro, list(lgb1$fao_cod, lgb1$choked), sum)
euros_choked_not_choked  <- euros_choked_not_choked/1e3 # convert in thousands euros
euros_choked_not_choked <- cbind.data.frame(euros_choked_not_choked,
                                            percentdiff= - euros_choked_not_choked[,"1"] /(euros_choked_not_choked[,"1"]+euros_choked_not_choked[,"0"])*100
)
table_choked <- cbind.data.frame(
  weight_choked_not_choked,
  euros_choked_not_choked[rownames(weight_choked_not_choked), ]
)
table_choked <- round(table_choked)
colnames(table_choked) <- c("tons before choke", "tons after choke", "% difference", "'000 euros before choke", "'000 euros after choke", "% difference")

table_csv <- cbind(species=rownames(table_choked), table_choked)
rownames(table_csv) <- 1:nrow(table_csv)
write.csv(table_csv, paste0(tables_plots_path, "/choke_table_France.csv"), row.names=F)

xx <- head(table_choked,20)
knitr::kable(xx)
tons before choke tons after choke % difference '000 euros before choke '000 euros after choke % difference
MNZ 8165 2199 -21 35788 12550 -26
HAD 3615 835 -19 7713 1997 -21
WHG 2641 471 -15 4583 1179 -20
HKE 1723 407 -19 4243 1364 -24
LEZ 1689 431 -20 5063 1760 -26
GUR 1472 302 -17 1039 216 -17
CTC 1313 485 -27 4352 1308 -23
BIB 1173 230 -16 799 225 -22
RJN 1047 275 -21 1957 513 -21
SYC 1024 194 -16 447 105 -19
MEG 997 216 -18 2917 875 -23
SDV 885 271 -23 1222 345 -22
JOD 846 175 -17 8994 1989 -18
COE 708 272 -28 681 221 -24
LEM 498 92 -16 2355 514 -18
RJM 463 104 -18 1038 252 -20
RJH 345 154 -31 787 345 -30
WIT 321 117 -27 706 277 -28
SQZ 309 233 -43 2739 1323 -33
COD 280 76 -21 1166 322 -22
# Table 3.1.f.3.

# looking at trips
lgbu <- lgb1[!duplicated(lgb1$trip_id),]
lgbu$counttrip <- 1:nrow(lgbu)
# X11(width=10, height=5)
par(mar=c(5,5,1,1))
plot(lgbu$counttrip, lgbu$choked, type="l", axes=FALSE, lwd=3, ylab="", xlab="# Trips")
axis(2, las=2, at=c(0,1), label=c("not choked","choked"))
axis(1)

plot(lgbu$counttrip, lgbu$cumsum_cod/1e3 , type="l", axes=FALSE, lwd=1, ylab="Accumulated cod catches (t)", xlab="Number of trips")
abline(h=294*fra_cod_dif$per, lty=3)
axis(2, las=2)
axis(1)
abline(v=lgbu[lgbu$month==1, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==2, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==3, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==4, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==5, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==6, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==7, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==8, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==9, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==10, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==11, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==12, "counttrip"][1], lty=2)
box()

plot_fin <- recordPlot()
png(paste0(tables_plots_path, "/choke_plot_France.png"), width=6500, height=5300, res=1200)
plot_fin
dev.off()
## png 
##   2

Economic loss if choke in 2019 for Ireland

irl_cod <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="IRL" & fao_cod=="COD" & year==2019), sum)
irl_cod_sub <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="IRL" & fao_cod=="COD" & year == 2019 & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod), sum)
names(irl_cod) <- c("year", "catch_all_metiers")
names(irl_cod_sub) <- c("year", "all_OTB_OTT_PTB")
irl_cod_dif <- merge(irl_cod, irl_cod_sub)
irl_cod_dif$per <- irl_cod_dif$all_OTB_OTT_PTB/irl_cod_dif$catch_all_metiers
knitr::kable(irl_cod_dif)
year catch_all_metiers all_OTB_OTT_PTB per
2019 580537.6 367775.6 0.6335087
lgb1 <- data [data$year=="2019" & data$MS_cod=="IRL" & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod,]
lgb1 <- orderBy(~month, data=lgb1)
lgb1$is_cod     <- 0
lgb1$cumsum_cod <- 0
lgb1[lgb1$fao_cod=="COD", "is_cod"] <- 1
lgb1 [lgb1$is_cod==1, "cumsum_cod"] <- lgb1 [lgb1$is_cod==1, "catch_volume_kg"]
lgb1 [, "cumsum_cod"] <- cumsum(lgb1$cumsum_cod)
lgb1 [, "choked"] <- 0
lgb1 [lgb1 [, "cumsum_cod"] > 461000*irl_cod_dif$per, "choked"] <- 1

weight_choked_not_choked <- tapply(lgb1$catch_volume_kg, list(lgb1$fao_cod, lgb1$choked), sum)
weight_choked_not_choked <- weight_choked_not_choked/1e3 # convert in tons
weight_choked_not_choked <- orderBy(~ -0, weight_choked_not_choked)
weight_choked_not_choked <- cbind.data.frame(weight_choked_not_choked,
                                             percentdiff= - weight_choked_not_choked[,"1"] /(weight_choked_not_choked[,"1"]+weight_choked_not_choked[,"0"])*100
)
euros_choked_not_choked  <- tapply(lgb1$catch_value_euro, list(lgb1$fao_cod, lgb1$choked), sum)
euros_choked_not_choked  <- euros_choked_not_choked/1e3 # convert in thousands euros
euros_choked_not_choked <- cbind.data.frame(euros_choked_not_choked,
                                            percentdiff= - euros_choked_not_choked[,"1"] /(euros_choked_not_choked[,"1"]+euros_choked_not_choked[,"0"])*100
)
table_choked <- cbind.data.frame(
  weight_choked_not_choked,
  euros_choked_not_choked[rownames(weight_choked_not_choked), ]
)
table_choked <- round(table_choked)
colnames(table_choked) <- c("tons before choke", "tons after choke", "% difference", "'000 euros before choke", "'000 euros after choke", "% difference")

table_csv <- cbind(species=rownames(table_choked), table_choked)
rownames(table_csv) <- 1:nrow(table_csv)
write.csv(table_csv, paste0(tables_plots_path, "/choke_table_Ireland.csv"), row.names=F)

xx <- head(table_choked,20)
knitr::kable(xx)
tons before choke tons after choke % difference '000 euros before choke '000 euros after choke % difference
NEP 4624 596 -11 27560 3453 -11
ANF 1367 640 -32 4949 2511 -34
WHG 1169 396 -25 1805 658 -27
HAD 1100 504 -31 1990 1022 -34
LEZ 1055 298 -22 3099 962 -24
HKE 519 135 -21 1444 467 -24
RAJ 355 167 -32 557 284 -34
COD 292 76 -21 922 256 -22
WIT 222 79 -26 387 142 -27
LEM 149 33 -18 438 103 -19
JOD 138 21 -13 772 135 -15
SQE 111 96 -46 301 297 -50
PLE 107 42 -28 218 93 -30
LIN 81 22 -22 128 39 -23
SOL 80 39 -33 779 436 -36
POL 75 12 -14 184 31 -14
CTC 74 286 -80 368 1327 -78
TUR 46 13 -23 456 142 -24
SYC 41 30 -43 23 19 -45
OTH 29 NA NA 70 NA NA
# Table 3.1.f.3.

# looking at trips
lgbu <- lgb1[!duplicated(lgb1$trip_id),]
lgbu$counttrip <- 1:nrow(lgbu)
# X11(width=10, height=5)
par(mar=c(5,5,1,1))
plot(lgbu$counttrip, lgbu$choked, type="l", axes=FALSE, lwd=3, ylab="", xlab="# Trips")
axis(2, las=2, at=c(0,1), label=c("not choked","choked"))
axis(1)

plot(lgbu$counttrip, lgbu$cumsum_cod/1e3 , type="l", axes=FALSE, lwd=1, ylab="Accumulated cod catches (t)", xlab="Number of trips")
abline(h=461*irl_cod_dif$per, lty=3)
axis(2, las=2)
axis(1)
abline(v=lgbu[lgbu$month==1, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==2, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==3, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==4, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==5, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==6, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==7, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==8, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==9, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==10, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==11, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==12, "counttrip"][1], lty=2)
box()

plot_fin <- recordPlot()
png(paste0(tables_plots_path, "/choke_plot_Ireland.png"), width=6500, height=5300, res=1200)
plot_fin
dev.off()
## png 
##   2

Economic loss if choke in 2019 for Belgium

bel_cod <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="BEL" & fao_cod=="COD" & year==2019), sum)
bel_cod_sub <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="BEL" & fao_cod=="COD" & year == 2019 & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod), sum)
names(bel_cod) <- c("year", "catch_all_metiers")
names(bel_cod_sub) <- c("year", "all_OTB_OTT_PTB")
bel_cod_dif <- merge(bel_cod, bel_cod_sub)
bel_cod_dif$per <- bel_cod_dif$all_OTB_OTT_PTB/bel_cod_dif$catch_all_metiers
knitr::kable(bel_cod_dif)
year catch_all_metiers all_OTB_OTT_PTB per
2019 29953.4 6316.052 0.2108626
lgb1 <- data [data$year=="2019" & data$MS_cod=="BEL" & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod,]
lgb1 <- orderBy(~month, data=lgb1)
lgb1$is_cod     <- 0
lgb1$cumsum_cod <- 0
lgb1[lgb1$fao_cod=="COD", "is_cod"] <- 1
lgb1 [lgb1$is_cod==1, "cumsum_cod"] <- lgb1 [lgb1$is_cod==1, "catch_volume_kg"]
lgb1 [, "cumsum_cod"] <- cumsum(lgb1$cumsum_cod)
lgb1 [, "choked"] <- 0
lgb1 [lgb1 [, "cumsum_cod"] > 18000*bel_cod_dif$per, "choked"] <- 1

weight_choked_not_choked <- tapply(lgb1$catch_volume_kg, list(lgb1$fao_cod, lgb1$choked), sum)
weight_choked_not_choked <- weight_choked_not_choked/1e3 # convert in tons
weight_choked_not_choked <- orderBy(~ -0, weight_choked_not_choked)
weight_choked_not_choked <- cbind.data.frame(weight_choked_not_choked,
                                             percentdiff= - weight_choked_not_choked[,"1"] /(weight_choked_not_choked[,"1"]+weight_choked_not_choked[,"0"])*100
)
euros_choked_not_choked  <- tapply(lgb1$catch_value_euro, list(lgb1$fao_cod, lgb1$choked), sum)
euros_choked_not_choked  <- euros_choked_not_choked/1e3 # convert in thousands euros
euros_choked_not_choked <- cbind.data.frame(euros_choked_not_choked,
                                            percentdiff= - euros_choked_not_choked[,"1"] /(euros_choked_not_choked[,"1"]+euros_choked_not_choked[,"0"])*100
)
table_choked <- cbind.data.frame(
  weight_choked_not_choked,
  euros_choked_not_choked[rownames(weight_choked_not_choked), ]
)
table_choked <- round(table_choked)
colnames(table_choked) <- c("tons before choke", "tons after choke", "% difference", "'000 euros before choke", "'000 euros after choke", "% difference")

table_csv <- cbind(species=rownames(table_choked), table_choked)
rownames(table_csv) <- 1:nrow(table_csv)
write.csv(table_csv, paste0(tables_plots_path, "/choke_table_Belgium.csv"), row.names=F)

xx <- head(table_choked,20)
knitr::kable(xx)
tons before choke tons after choke % difference '000 euros before choke '000 euros after choke % difference
SOL 28 36 -56 350 446 -56
RJH 27 31 -53 55 63 -53
RJC 24 19 -44 41 33 -44
PLE 17 20 -55 41 51 -55
LEZ 10 18 -63 20 34 -63
RJI 10 4 -29 16 7 -29
HAD 8 2 -17 14 3 -17
RJN 5 0 -3 7 0 -3
ANF 5 11 -67 51 104 -67
SYC 4 8 -67 2 4 -67
COD 4 3 -40 12 8 -40
WHG 3 14 -84 3 15 -84
WIT 2 4 -61 4 6 -61
POL 2 1 -24 9 3 -24
DAB 2 3 -65 1 3 -65
JOD 2 1 -47 11 10 -47
NEP 1 1 -34 7 4 -34
GUU 1 NA NA 1 NA NA
LEM 1 1 -62 3 5 -62
SYT 1 1 -49 0 0 -49
# Table 3.1.f.3.

# looking at trips
lgbu <- lgb1[!duplicated(lgb1$trip_id),]
lgbu$counttrip <- 1:nrow(lgbu)
# X11(width=10, height=5)
par(mar=c(5,5,1,1))
plot(lgbu$counttrip, lgbu$choked, type="l", axes=FALSE, lwd=3, ylab="", xlab="# Trips")
axis(2, las=2, at=c(0,1), label=c("not choked","choked"))
axis(1)

plot(lgbu$counttrip, lgbu$cumsum_cod/1e3 , type="l", axes=FALSE, lwd=1, ylab="Accumulated cod catches (t)", xlab="Number of trips")
abline(h=18*bel_cod_dif$per, lty=3)
axis(2, las=2)
axis(1)
abline(v=lgbu[lgbu$month==1, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==2, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==3, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==4, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==5, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==6, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==7, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==8, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==9, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==10, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==11, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==12, "counttrip"][1], lty=2)
box()

plot_fin <- recordPlot()
png(paste0(tables_plots_path, "/choke_plot_Belgium.png"), width=6500, height=5300, res=1200)
plot_fin
dev.off()
## png 
##   2

Economic loss if choke in 2019 for Spain

esp_cod <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="ESP" & fao_cod=="COD" & year==2019), sum)
esp_cod_sub <- aggregate(catch_volume_kg~year, subset(data, MS_cod=="ESP" & fao_cod=="COD" & year == 2019 & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod), sum)
names(esp_cod) <- c("year", "catch_all_metiers")
names(esp_cod_sub) <- c("year", "all_OTB_OTT_PTB")
esp_cod_dif <- merge(esp_cod, esp_cod_sub)
esp_cod_dif$per <- esp_cod_dif$all_OTB_OTT_PTB/esp_cod_dif$catch_all_metiers
knitr::kable(esp_cod_dif)
year catch_all_metiers all_OTB_OTT_PTB per
2019 15066.54 14165.64 0.9402052
lgb1 <- data [data$year=="2019" & data$MS_cod=="ESP" & data$metier_dcf_6_cod %in% all_OTB_OTT_PTB_mets_with_cod,]
lgb1 <- orderBy(~month, data=lgb1)
lgb1$is_cod     <- 0
lgb1$cumsum_cod <- 0
lgb1[lgb1$fao_cod=="COD", "is_cod"] <- 1
lgb1 [lgb1$is_cod==1, "cumsum_cod"] <- lgb1 [lgb1$is_cod==1, "catch_volume_kg"]
lgb1 [, "cumsum_cod"] <- cumsum(lgb1$cumsum_cod)
lgb1 [, "choked"] <- 0
lgb1 [lgb1 [, "cumsum_cod"] > 0*esp_cod_dif$per, "choked"] <- 1

weight_choked_not_choked <- tapply(lgb1$catch_volume_kg, list(lgb1$fao_cod, lgb1$choked), sum)
weight_choked_not_choked <- weight_choked_not_choked/1e3 # convert in tons
weight_choked_not_choked <- orderBy(~ -0, weight_choked_not_choked)
weight_choked_not_choked <- cbind.data.frame(weight_choked_not_choked,
                                             percentdiff= - weight_choked_not_choked[,"1"] /(weight_choked_not_choked[,"1"]+weight_choked_not_choked[,"0"])*100
)
euros_choked_not_choked  <- tapply(lgb1$catch_value_euro, list(lgb1$fao_cod, lgb1$choked), sum)
euros_choked_not_choked  <- euros_choked_not_choked/1e3 # convert in thousands euros
euros_choked_not_choked <- cbind.data.frame(euros_choked_not_choked,
                                            percentdiff= - euros_choked_not_choked[,"1"] /(euros_choked_not_choked[,"1"]+euros_choked_not_choked[,"0"])*100
)
table_choked <- cbind.data.frame(
  weight_choked_not_choked,
  euros_choked_not_choked[rownames(weight_choked_not_choked), ]
)
table_choked <- round(table_choked)
colnames(table_choked) <- c("tons before choke", "tons after choke", "% difference", "'000 euros before choke", "'000 euros after choke", "% difference")

table_csv <- cbind(species=rownames(table_choked), table_choked)
rownames(table_csv) <- 1:nrow(table_csv)
write.csv(table_csv, paste0(tables_plots_path, "/choke_table_Spain.csv"), row.names=F)

xx <- head(table_choked,20)
knitr::kable(xx)
tons before choke tons after choke % difference '000 euros before choke '000 euros after choke % difference
LEZ 5 2274 -100 19 8586 -100
ANF 3 2813 -100 NA NA NA
RJO 2 115 -98 NA NA NA
RJN 1 130 -99 2 327 -99
HKE 0 1890 -100 1 7471 -100
HAD 0 81 -100 0 124 -100
EOI 0 90 -100 0 237 -100
LEM 0 67 -100 0 179 -100
A NA NA NA NA NA NA
ALB NA NA NA NA NA NA
ALF NA 0 NA NA NA NA
ALV NA NA NA NA NA NA
ANE NA NA NA NA NA NA
ANG NA NA NA NA NA NA
API NA NA NA NA NA NA
ARG NA NA NA NA NA NA
ARU NA NA NA NA NA NA
ASD NA NA NA NA NA NA
ATB NA NA NA NA NA NA
BAR NA NA NA NA NA NA
# Table 3.1.f.3.

# looking at trips
lgbu <- lgb1[!duplicated(lgb1$trip_id),]
lgbu$counttrip <- 1:nrow(lgbu)
# X11(width=10, height=5)
par(mar=c(5,5,1,1))
plot(lgbu$counttrip, lgbu$choked, type="l", axes=FALSE, lwd=3, ylab="", xlab="# Trips")
axis(2, las=2, at=c(0,1), label=c("not choked","choked"))
axis(1)

plot(lgbu$counttrip, lgbu$cumsum_cod/1e3 , type="l", axes=FALSE, lwd=1, ylab="Accumulated cod catches (t)", xlab="Number of trips")
abline(h=0*esp_cod_dif$per, lty=3)
axis(2, las=2)
axis(1)
abline(v=lgbu[lgbu$month==1, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==2, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==3, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==4, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==5, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==6, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==7, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==8, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==9, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==10, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==11, "counttrip"][1], lty=2)
abline(v=lgbu[lgbu$month==12, "counttrip"][1], lty=2)
box()

plot_fin <- recordPlot()
png(paste0(tables_plots_path, "/choke_plot_Spain.png"), width=6500, height=5300, res=1200)
plot_fin
dev.off()
## png 
##   2