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)
| 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)
| 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)
| 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)
| 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)
| 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)
| 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)
| 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)
| 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