From 9634ec12f9b0e18ca8447da96212ebc5019317a0 Mon Sep 17 00:00:00 2001
From: Benjamin Blanz <rganb@gmx.de>
Date: Tue, 20 Aug 2024 15:25:33 +0200
Subject: [PATCH] Add documentation and GTAP capability

---
 .gitignore                                    |   1 +
 scenarioProcessing/a1ProcessScenarioData.R    |  10 +++++++
 scenarioProcessing/aggregateShocks.R          |   4 ++-
 scenarioProcessing/funAggregateNuts2CNT.R     |  19 +++++++++++--
 scenarioProcessing/funRelData.R               |  26 +++++++++++++++++-
 .../helperData/GTAP_NACE_sector_mapping.xlsx  | Bin 0 -> 12966 bytes
 scenarioProcessing/relativizeShocks.R         |  17 +++++++++---
 scenarioProcessing/stockAggregation.R         |  15 ++++++----
 8 files changed, 79 insertions(+), 13 deletions(-)
 create mode 100644 scenarioProcessing/a1ProcessScenarioData.R
 create mode 100644 scenarioProcessing/helperData/GTAP_NACE_sector_mapping.xlsx

diff --git a/.gitignore b/.gitignore
index 82f9275..ab5c6e5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -160,3 +160,4 @@ cython_debug/
 #  and can be added to the global gitignore or merged into this file.  For a more nuclear
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
 #.idea/
+.~lock.*#
diff --git a/scenarioProcessing/a1ProcessScenarioData.R b/scenarioProcessing/a1ProcessScenarioData.R
new file mode 100644
index 0000000..928275c
--- /dev/null
+++ b/scenarioProcessing/a1ProcessScenarioData.R
@@ -0,0 +1,10 @@
+#
+# script to perform all the processing on the scenarios we got from James
+#
+# Benjamin Blanz 2024
+# 
+
+source('stockAggregation.R')
+source('sectorMappingNACE2GTAP.R')
+source('aggregateShocks.R')
+source('relativizeShocks.R')
diff --git a/scenarioProcessing/aggregateShocks.R b/scenarioProcessing/aggregateShocks.R
index bde754e..95b273d 100644
--- a/scenarioProcessing/aggregateShocks.R
+++ b/scenarioProcessing/aggregateShocks.R
@@ -1,6 +1,8 @@
 #
 # aggregate the flood2010 and earthquake scenarios to country level
-# 
+#
+# Benjamin Blanz 2024
+#
 
 source('funAggregateNuts2CNT.R')
 codes <- read.csv("helperData/nuts3fid4Codes.csv")
diff --git a/scenarioProcessing/funAggregateNuts2CNT.R b/scenarioProcessing/funAggregateNuts2CNT.R
index 90028f3..da9f1ee 100644
--- a/scenarioProcessing/funAggregateNuts2CNT.R
+++ b/scenarioProcessing/funAggregateNuts2CNT.R
@@ -1,13 +1,28 @@
+#
+# Function to aggregate stocks or stock data to the country level from the NUTS3 level
+# 
+# Parameters
+# 	data		data to aggregate
+# 	codes   country codes
+#
+# Returns
+#   aggregated data
+#   
 aggregateNUTS3ToCountry <- function(data,codes){
-	sectorCols <- grep('ALL|TOTAL|^[A-Z]$',names(data),perl = T)
+	sectorColPattern <- 'ALL|TOTAL|^[A-Z]$|AGR|MIN|MFG|EGW|CNS|TRD|OTP|WTP|CMN|OFI|OBS|REA|PUB|OSG'
+	# identify columns with data rather than identifiers
+	sectorCols <- grep(sectorColPattern,names(data),perl = T)
+	# convert data to numeric (deals with in import error)
 	for(s in sectorCols){
 		suppressWarnings(data[[s]] <- as.numeric(data[[s]]))
 	}
 	data[is.na(data)] <- 0
+	# combine the data and the codes datasets on the fid4 field if the data dataset 
+	# does not already have a CNTR_CODE field that can be used for aggregation
 	if(!('CNTR_CODE'%in%names(data))){
 		data <- merge(codes,data,by = 'fid4')
 	}
-	sectorCols <- grep('ALL|TOTAL|^[A-Z]$',names(data),perl = T)
+	sectorCols <- grep(sectorColPattern,names(data),perl = T)
 	data <- aggregate(data[,c(sectorCols)],
 											 by = list(Category=data$CNTR_CODE),FUN=sum)
 	names(data)[1] <- 'CNTR_CODE'
diff --git a/scenarioProcessing/funRelData.R b/scenarioProcessing/funRelData.R
index 1b3cb75..981c275 100644
--- a/scenarioProcessing/funRelData.R
+++ b/scenarioProcessing/funRelData.R
@@ -1,11 +1,34 @@
+# 
+# Function to calculate the relative magnitude of shocks to capital stocks
+# relative to the size of the sectoral capital stocks.
+# 
+# Stocks can be provided either at the country or at the NUTS3 level. The regional 
+# aggregation of data and stocks has to match. 
+# The column names in data and stocks should match.
+# 
+# Parameters
+# 	data  				The shock impacts specifying the capital destroyed
+# 	stocksNUTS3   The sectoral stocks of capital at NUTS3 level
+#   stocksCNT     The sectoral stocks of capital at country level
+# 
+# Returns
+# 	shock data relative to the stocks
+# 
+# Benjamin Blanz 2024
+#
 relData <- function(data,stocksNUTS3=NULL,stocksCNT=NULL){
+	sectorColPattern <- 'ALL|TOTAL|^[A-Z]$|AGR|MIN|MFG|EGW|CNS|TRD|OTP|WTP|CMN|OFI|OBS|REA|PUB|OSG'
+	# Ensure the Total clumn is called TOTAL, not ALL as in some scenarios
 	names(data)[names(data)=='ALL'] <- 'TOTAL'
 	names(stocksNUTS3)[names(stocksNUTS3)=='ALL'] <- 'TOTAL'
 	names(stocksCNT)[names(stocksCNT)=='ALL'] <- 'TOTAL'
-	sectorCols <- grep('ALL|TOTAL|^[A-Z]$',names(data),perl = T)
+	# Identify the columns with data (not the columns with country or idx). 
+	sectorCols <- grep(sectorColPattern,names(data),perl = T)
+	# prepare an empty data set for the relative data, keeping the index columns
 	data.rel <- data
 	data.rel[,sectorCols] <- NA
 	for(i in 1:nrow(data)){
+		# identify the correct row in the stocks 
 		if('fid4' %in% names(data)){
 			stocks <- stocksNUTS3
 			rowInStocks <- which(stocks$fid4 == data$fid4[i])
@@ -22,6 +45,7 @@ relData <- function(data,stocksNUTS3=NULL,stocksCNT=NULL){
 		if(length(rowInStocks)==0){
 			warning(paste('no region match for data row',i))
 		}
+		# calculate relative impact for all sectors of this row
 		for(s in names(data)[sectorCols]){
 			suppressWarnings(data.rel[i,s] <- as.numeric(data[i,s]) / as.numeric(stocks[rowInStocks,s]))
 			if (is.nan(data.rel[i,s])|is.na(data.rel[i,s])){
diff --git a/scenarioProcessing/helperData/GTAP_NACE_sector_mapping.xlsx b/scenarioProcessing/helperData/GTAP_NACE_sector_mapping.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..d2160144f58b37fd0c309d8698b79a23b824b104
GIT binary patch
literal 12966
zcmWIWW@Zs#;Nak3h{_U)WIzH^3=9kvIr{NMsX4{^<@rU~N%{HNdKI}jdLWUtL5^96
z4FvW)7fpJ<_05u&Af4wzTsxP&WNz8B%_qk*!#wxX?^V7TG8_xp=N~zfmEM+k`fA~e
zqv~$2_5=r{SaOscyFNEtwC2sdd2tD$XHV&Vn{`=yQ>?<oG8xScXC8)n8Hus`&g6e_
zY01aJ3v&{3QZJ;fkDl_`CjOMgwT4%2!8uv1g~8$qWy>4>U;O5KSu&=&SmMG*S?-gW
z;*%Hj>27L@uTbmO>1^CHKhZ9g?NM-l32))C{x{;mx9WfHQ13o*ouxj6>-%S$H;*bW
z>8<ID*k9~(+;m-vPyC@XUX{f~SIq;w**W%Iu#$Yhz`(GQiGcxkIPpWm2|0v%gHPt&
zHW1kRU7O`zuIIfh>1oX2+a@r?Hd*+dn_{vlB6x{I-*%C|-+MH)eJ{#YZ8#bE<Js9W
zmmfWg@qEv5jBBlu*f-{9ji*ZVuIMG*m!EuZN~PG&Q#-G4DU~sv&YLDLrk?w4<+bYm
zd)m?)g*MLnz3;5oiXCrXONqF0<!be%%$l2@&1|z>l_BhsJj=l&>vsgKah);y|G^8c
zibfHy-vl%s5<AlH#N$n_X4dIx-A*1AS<;^>V(XuZUhLa>X2$-S#kKi+rlkLNo*Phh
zgiB91VDVDBTJL&?wkv|uCskUybEF?!;r)}@NjWu0vdFAYol*U_^}QK<AH$?SE)p%;
zdA4VoQBnrq^(+?=*Kg}jY?HjA?4SGG#-4lX>)<09uG7oEPu^|UE#n>aAjLwqbYHZ<
zZXsne-gM??ivIGm99?r{E^Ky^=w8?1{P4<2!4+}IfzJPS_uV<S>vC}3f_JGm0+xQ(
zTYkz<-DBb64N6mw+&II^dM4Jcq(hcx>+FkvSeq99iI5iw<dF7b>(aMU?+{%W`#Asf
z&Rf|NPIRW&b2M!@*lJw2?Zvs(Q<g0I9bYs_$@Bu(iOkMQc4eh+pC6|fbS|Iln7@Im
zeXqftp6dB+iwpH{s&M-1#U0=LMEvEch38(1Z~vn!^~*zU!Hu|#?aw=M&q<z``uIt3
z)bIBze$G{!wjPvj_HbmZ`OV0{P{52Y-3UX{O-V*-ZmK?*Gz8_Yqv5yn<`@Xnz7LPs
zW6n11wbtwvS-W02>jzD)Yq&ncJDF$G+$F!i%kn(jwBe=K-6aYeD{7vd*_fWQdHVJJ
z-hFa2LtItOrWGefcW_3jww^G%c7FHU-zS{BgeMBbac<lA%_s7D_<OyRQP0<;aNKE)
zn%p}zJ9fQcn8+EI*br0U3)90ksYuN-I^3z^?+{(|=2lsQcM!+w&(oLx2+xcPC=0&W
zK8f@0B#&dAY%jI+vr-iuy)OK?%WyGP`NFQMX}>x*AM-b5JM~o8^w)7Yw%K>A&-PEc
z`A&KIi-{k4Q)iz(6h3_idx7oorO#d!$SUkCe_nH@*Jk<i;^=*reXeKkewJidDB@<=
z|Mo|t>`vJ-maKU39gi<W`hHaYcXf_EU)ehIgX`J$n9C<iZ9X#p!CuDq=RX`NEB~0)
z^qsxN;Ba+Mi{9)9yc*A2cg$bhY|z!jch~t+SAAT|amR#qPj>l$iV(}S`!5`SwL)dT
z)=t(2Pd$;94ps4+rvx0RT~WI9g73mCzi(ozI#-Cy$h}^;;DXKVxB2{4+`gZ6HM#e(
zceG~sElH9Jy2PQ9?dEo7ZRaXxVLO2ku}Ar@_6cgPUcOz({Kdh?m8P3x`lMN(Z?8Ik
zQ#^Ic>Mrq<Vm-?e6Kcg?%f#s(zVUs_+Ks_#CqJySy5+NQSJ2sxxCgB>t2OqAd2V-@
z^W2xoY>v~5Yf>$T&l!qNTv@1a-(tSpnQ!7vt*b9LS$?xz<a*Bg;)e{g#IFjMWHk0&
zcXR6!IyWh;&N}E^oxv(WbJrzjgJ(To$NaF_XRda3z_HH@kKWTx&D9W`I=Mx5*7-Tx
zE#kPu-Tv$5SL|!pIb+)e;hxL9vlFIVDOh&s9n1Pu-oU9RbX+Rf;<ZAyA6~oklDW&%
zw-?SF%3C_MHBIQ_q82CF^ZMe8Eeu7zNgpYk{j>V0#g781$XQNdeGzSHU$&}r?|shx
zST%KbPW3n2zaniXtG@T^c5apxU7WCbX{e<8-qOnYU*Fees3=BFNsymda4b1hM6Y7v
zx2@k79^Q0#GC!})@5>)@cK%`p<)9D8RUR@jGcXwQ;>$t2kQ`K8Qkj!l49YuiBVzMs
z8;I22U%%ncKBu5r>x37+vP!+>9lCa<N#?C)X0|m_(l`C2f4wg6ZFKmx(X2XR)ichT
zn)7?xANO1E`G>Q*Zx6{-J1Wx8SJShiBlVDr+0pBtZyxuT$aHpb+96WXrM>E7<`=~)
z%N~8vT7E_0gu)k*=f|eZxAi}6soP;UK~?m^WUo*GX5XomyF#X@+?vGouD7=Sk!TC2
zOlwl)8Lq{FI@1p|^zb>o`)JEpAe8Xr(%jQ5OKKK3X&u@q-WoSqfh$<^T4mn`=0ySb
z7e^mbJb9w+_$2RE9*(8kBovxNuDH&b>XudFIAv1#Hg#3j$5UTRUw<-3gLnShPgk$T
z_V#-HbBjAV=USap{qYyv*AuTyy0vfm?lM1%XpL`HhhMeJ8Xv1PPfLE&prO-g7}6<e
z67bn>z1Ire1hrizrFY^N8N05!_+y3p3i*I?#ZV^W{r0bAT7=w;G&%O`Ew@mO3BKbP
z|8#D~L5(LSrV}DI>TLH{_CBv>=qC2zQc+Ik^qtj>pA4V=Ii<IF{>xvlYdv;($@^6Y
znkSS-y*<17^c542$WK{nxs$4X$b4AJC+UzrCw0*R>Bi|^Qc|<4#8*7jP{~y|wKVDb
zgj;+G*Vb6=&Xf>~K2fXh*0$~pe`d#NrdjVqG@mn+sTobY(z?+9c~|G`opUy>lg`ny
zmh3qFqfKq?i8-ncZ(XZammJote|tG2dcM@$*=}DS?hAfvQemBRtlH4O<7QRv+oNxw
zFh^O<;?C0eQwdvwvzq1x8Q)ViRsSE^bKQR5?!z_jW#n_Jd*YY!PyPP&&*CrJ%a+OC
z5ZiNZ+jRaP!vCA+nd`DOtlF2Z`>7{=-KFZSZzfi5Yo52asC1w5`8{dVU*C-_vA)Oh
zJs>5&ZP|}Q_pW`Mm)ouSxQ+ArCpn>Wseu#B>Pk!$Wb8a0_W!vrn`xxrE2pYEb&=9K
zt8*V4gl&#XUioG6K-T`)?j@)A!rd(;1MaN7a?2-Ujegi-dr;muIJx)wLna1>xuW>;
zh6E&UfLbNR8L6oy#rj|hT*d5-jLm;!AyRukUhrux$2su;#>IRW=1#S`eR11n<+rnE
zOKg~RB=We!?fCuI{XCj;XJolnW(xm)Qa7hMCpPxj{P~~H+!Trn4D*S(@_IL)rqHxY
zNv4JMHg{^IcKN9mOnlTMn(}Dhgj0WQ{^ssU{<AVlSV_hpAVj!a*QN4lW#}d~%XL+2
zR9N5WNS@!=5Werzs!76gR`q{6^ZEDR8>bC*@7y{hTHnvRE#^|onaE8iJr%Q*=cGh9
z{J1)6+SPlviht=<YaX}uD3$kDKT^5(@-63uiTcYq(?z4Z;tqXMGci6r(cJ39=b69%
zpBL&<xT-sK>Zk9=K0UFGd^W>A_mqFaZkA<onQlkBy>86yoBH?toybc*0{?oK20T$c
zr^QsQ8aJ_8scuSEyv>I98#dhkwk&!5&k8F)%`47s&k8op*>U~fRi>jtiv;R}H?6ih
z@3%E=w=_q=<cC@L{<{CpWmb7}Gp%mQpLKNqjl5%<);qAh*s}X0KjXD;+;i4l6L_Q>
z*>!Enf&T9*%AE?C-*@lsJ#paVU+;e>HeGY@npQEHQ<7;%a?SpItsUxlia*aBP|;r6
zEqT)Dgz!#3rH%Cse=SPph9CR%n<f2#(5vs(@9S5c{~CJoxo%yew5#`K?fz@Z$^J#2
zx89!CTB)P5Uv=8?4V{Nd4c6p&%AQqvx9|FVtE>B0DRuu?n{cV%du9AAX;sM;rH5HI
ze3-l}$+ob4uLSRf0)eeMj=uYZPFk?=u3jI|`M3J8?P||`s?n1eqgXG;GQTd(yVv;n
zr3A;?EAgJWZKu7&)*H^s(K`4^ZPhND>{-i}Y}gfQSNw;^@KAl!ivLq9CO)0$ay{)&
z-FbP>vU9I5<;VOgo6cNPH{J2E(S{u|-`Xzzv~BA>e&LbfABnid$5K-so4!~Wz`0i}
zuy%5C<MQNlWqjdRSIAzk%C7Z%_-o@6kN!n%^98nZ9XP>%Y1+SasXwQ_d+?YmZ2RSZ
zH*TNekUk<V+{h<rEc{Sj>6x06g^P}l%MB$Ph2^Ol76%2atZf^)PQPyv&HQquV?lF=
zX_e0rhsDQTctDI3e*P7ijIC_W9^nBQEd3>%kqW!obOd5u#7+o`x{8MerW7qXtzqPR
zy)>oqh10LCQ(IY$wiJuRM|o}ZyzH@9(`?aq)2(N91blOgxm|*9A7bKsW5TPawY2fz
zMENWk-ZK{<vg)7Xqa$ZbNDh*Ui}G@wtoC+IjpQuOtnKg4FkNJ^G*@g?JJHJ8^v0xD
zPiyJKgMka?Ht8!Hm%KeUH&S7i!gH|>JymX%ZN)9W)IGS|_P%&y<JoT&xNC;g8AqmF
zvvMS2{lt7$ygr%sykmjDt}h2C%b9L+S~9uyj7sfX9d5M?CpRTH&D9Yup0XzR?DwDj
z@}GQ{eOb5QpiIh&?W|@!_rz})RCQ+_lu3ORBK*%V{X)(W=e_0*X_IAMo^9k^T+VLR
zb5Ab9!2P7v)k*aYYtFRUaFr}O<7UZKx=d$=NBXIzi@$xjpJ@ay-}R|~=h`#pY>d{G
zZBFj%S$((3a9x?OuypwDPD`(AE5H3oKDXfJ{Nu;2=6yT6v2bg&!Lh4(2PMLHcN(7f
zWTiZ9(q*SRQ?tb0ui>0#ue7JPJUoAKZ0*x|9c3&>ENW9*58Uh#Y5H*PK!*2@1)Vz|
zUTuDrKV!Xy#kngCJMK%HC_j7_xgp5b@v7X+2BWw4uVt;%xv$jH#_{#9!neP59QSo3
zg}>>pm+i~<T7E?J;rB07SFmk<)l+}<!XC-Oh0nIi&fe8$Qh0po_bG3?+D$dzDR^Jt
zSw3OOuh8U2(|Wl3RyG_ko?{a${_?!3=v$e4o|f-=_6Sc(HJ)~QzPpy?Ry{tCxzk0s
zw485m{<&?^hDk!E%1v*w!tz8I{k*1oC+rb0pAw$G?$*=YCb2p9&mAl|b7_m*^=1dF
z^B*Fj*G#|te4TMcZhiEuf9#+FzWMIr<Fi>A80JgiE8s;S1$=QvVo_>Ja7j^SUOKo9
zG&ejp|FVHd?ep-6N~X+6wo44Txy!eGNp`rt?&XWkg?ZB4$N9Deht=6$=X9Ria>O+M
zr-kM5=l5P$&%JjxdFc(tHIv=GPrm+KoReACNAc^v`8GMrUjA;eJYRZn>ctEGeGewh
zKG;34c;3V1x4urk^!u0A={Xt;AHBPCWx>a@=WcXW@7enM%vG_9hh7ZvE8m<wdi3_1
zM?va7Gu(dHeR^k<*gNgMv_v`oy*Rbzd6S#tw2Pi*=><qf%G_X9X}+1ZBPlu9LH~1e
zoY;i_PR|wgY*hJkx4>rM=It*yzTLE{QCRWi@Xc9mx1%+so=ST^=}omw4sJfb*LI%a
z8g{d(>yBtCUi@IIYu0_+TFNhL)fbJX4PQ&|7i9N`R3%nqE%Ew(=VHq;W2wH^X~$XF
z=Ixqf$Qc^#e8@K~TFNi+X>zYef2!f*b*p=IzdM}%(f&W=r0s`Q!qQ&dyF<k<CG&Gf
zo|(NX%w5QEMKYJy1*V$mmkZn;>t(t`@V#u-<lGv0Ice+0Jq+K9qNP~-v@?>ou$^C8
z<(eI!a^?Py$knyI<)+7%%#e+i$mxE=JtatCSL?R^BZhY-x5#XqG;Ox1sokxs$EtE}
zeN$j9l>BT_B)94Qud^5R?q8l|<fi-g;p;2Pyt>>jmz@=TLT4_D4SVrq*$cDtcXmWo
ztt&`Bdi1eb+;(YcjlWBdwC(WZ+mYWn>oIGtO;4ckdErIDy_>#HOp7Z&K3hHbJWJ+P
ztEH?3TBg^XibT^N2dh*X6t0eKE)}0Y?P<p~3HRQ~QggaDZ}Q0WjE~&xT;#(LZXw^o
z^7g4;*Sd4&3a2J0lqS8K#oNurZeNz|b7$e}l*EcTCha#O)~Ii3zA9CqDQpv29hLq=
z=CsuA4_BGrF4>Z`PfPPr&<>yQm$f{<Gvl20&l2hVQG3&D?XpL+&R$%jzkA8dZC%@s
z)phILn`8EV+CCQP^y<0KUaHkDF!g%<F~s;{McyILRoXiDjOQ$7s-M&JKKq-`|4PCC
zJNbmm^Q71<W_Rb^@{@BGI2hHLtZ^hG^Lxzc`MTNck`M3STe+C$|68{G450kf$GX&w
zkC}l%k`G^g;$~oAfQ)s4$G7%gvGqS}Ai(xu?|!9E_jY(LRG6}a+v!zXz{H5&+MfYQ
z^Hbv$|5UHx3he57DNwql_WQkLe!Xjp<tCXI>~jv@q!YV=rN`y!%w;L(Z@2%{jSk=Z
zh?C86;gRigzk0qf+?S#|aqr_hr^U8>(2O#8xWsKaYwObN9ip4Wc|X)_j41qW+#KpY
z^Q7C<Gr1whCcSNME#9<ooi*oojdCfTxe^~&E2(f_XPBBJ^TT@AKhcf;&Ki%G>;C*X
zZU6g|rV|31Sf;;7$cueZcBRq2K;+c%4R61T8eS+mb-Z!n0gZsdNgK@6cFhV>DW4bc
z|Fg^QSrRX+sxEoY_~z?*`LN!5juNrUYxLHhOrAB{J9}FEtjL4?brW6*gMx8>fk7W1
zBLhPP3%+0!WME)O$xjX_$}cF^PtGq&1@+5%g97~z83@#U57()czuPurXV0Vs2gN;n
zj%3uVzH8W(eSHh-k@$S;1&cl;Tx#F_dtY67{Q9Raf=w>6ZF#s>X3>`E3piA=we?;r
zR?YK&B;$H4DDvX6sRl0A1{p5vyq0lY%c{7tmMPbA#tBZ(U;|SZ@9faVjfvg@SH8$N
zTNPZ)a$M_O+AEqWzlCj%>PcR)giKf6qj3gdj)Jl(+F~kKqs-?z95paeaENNY9M@I9
zLu`wIxPHkA&ftRG!37hmSeSob{@`#YZ~NZd{fl?~4E|#^+paV$S?_PKp~Dgzo*)C>
zD^8-8jY;?4b?o`_{@rt7mrW|W`zG7?<v;g~Tf$s^&+E_ged2eHuT=}yTP*o?PQa%E
z*0#)bYi_H&JP^-1-8|%4_V+g|DSr&!X)Jsqek&k3J?|?w-_;8DV+!FRJK~a$|I+$s
zsk}p^{g!=B+Q~cBXSyBBl0Iu(IJ+?M@ZWls%<Cpg?lnj5U*B+Z@4*Q>ZWP@<u&;*a
z{6CGKY@k@XefQ0Saz+M*_jt!J1(0Jcv7i7HX>%v|W*;`-X?y=$^x&6S7QB%!x*l&b
zEYr$qU-Kr+COXsh5r59oCI9R4PT!iRF7Tq;{nF{M$&#CxTUN}?J-J82lP%0@>DzUx
z^Y2FQQ;2j)oGLM=b()OgB+>c=`Gb`lD-_c@6Oz5Zb<}xgPN_Cf5?h%YyuAD%hfn5>
z!0X?m4y~TNn6Z0{*6FrGEeWD4?!-4faf*@qus`7Glb8E;yzI04Xz}{M4@ZRpxBSHC
z0X%cf4$KewIcr1udv;w`HP;_e_3Rb)>lcK^hPZyZ8FjnxR-Nvw+~D+b;oDyh^Ls2Q
z4Luj<X1iif(ige%q<p=k&p(p`)vWpdT$R}S@>D$|DBS$Ik}ZERFfhzv#FyGZV|a+v
zURqp|p9>1Qr4tY49Wvl?c>lXg<GYZ$?1e?!ID)haXP9K>d|%<~tvV^h^M%9Zf3=R%
z&em@>$?D!df4*#+%}QperKNi>S}qel;PL8~$=-KGfA*|i#G?GXHSfs5jG2!&msCXF
zN$tPJbLm%{l?cbJrxox2ZcsRKwR&aC?3OKjGp|Lc)^Oyy``rGvRmZnBB&hljL)sEM
ztNaj?n#Uj1^9_IPSCS5ne{<*bI~$g`iM0pMuL#|jU-T!;uQJjvQEH>_DP^t_*`<Y6
z{R)f;*~eCwY^^+=u|veTXG=h0ao!ewW4#Nx!YdZX%d~|*&t+8S;AxogsjK*0P~65D
z_l|QK@Ohs4F=t(x=yJ)OUsIlBef8h}<JEVC*>N_W)@{u>Q(rvsD|y%cJpc2_w=<7|
zn>h#d9@i@}GB60@O_`t&P6h{TL~f3LW=U!;xC^^<lB4$#1A*4}|GHGaoA^#%cKcG<
zL{@)2j&%=A=c?Uym_0@1^}F-iSFfuQ5ZQC~;okRqC70csv%K)@q8fv9W*uspQJUQA
z%)eb&u>7xV++yGA;PY4C++9;+qP0e>oA>IS1&YUZDybY;pe8x#WJCGEmopk>cU->J
zpcFZa%V?!WwWR9LMGN@PSkCaj@Uq0qPp9x-!N$N(dv@2?Bqt|3PZ7Iv)%5WFOTAa~
zrWNk^!1YO?cFU>!xxM##8Ml_ce_)|0U4H3Q^Ddv<z0;qs?J>C1_qaN$dBy|QFR6X3
z!6L=iL?guZ9NzMNp@2+={^zAp22)Lwn-|$@x!kH^iRJM9EE95GJoEVFO$&Om&o*rI
ziT}7qAeGPPa6lSU@yvDNdrlazAGEsm)Fin`b9Z<hSNQq=+vB+=-IvmvG~dr<-+$-d
zT%hFbmc#kOhLM566>s2!1`$vLA6%|OVql_U_F)AHaPhF)w9{?!r6oF^+<SXB9$a^9
zM%W6eD}4uV)X!hId(GnIvhh!AYxbQ#`^|a&<not$HZ5bjsm0o9a;Wc=O!kLgo7z=x
znDNiPxZYiE2miH`&J(OUO@+_b@cJp9$!^(W;d%D&sdJZ(oIYuF*16|h#~Tw(#cFN6
zs96o4OqP1YD0%OHY3!jC?fPWl`h)8C%2Tdy<(VdVaKe*Z&edh{4tK<BXTCl4ZDMKs
zg2<Y6{C1~L{_o>{|LBUm;o6qX_xuu7*ZeEowJSsNhmq{8*e&Ibp8OZ=7o74bd-ib4
z1<9&MdJ-?KvYR$P|DEE)DV!YY(3Q%@DJSuyhsQwL{h)^6%@f_bQx9c?a?EoMd@3R|
zZ`wwmVh@(1%yL~i$7?hf&px{Fc7N_ZaB*_)*D6;nMg|5Uyy1<R*o{ENN$$Cgd`$)d
z$3EU&u5fpulz6gN<b;zJ51A7zmvFwFc(3!*H?hUdoS)x+`tK*IZ>}QDvgt+f#?#xB
z81px5{Wq^WJ9a61;EP@@nZ&FM#ue+H3r=!c&DFg6V6XFG2hTSZr(Dn5$Sqa5+?=gZ
ze2qorjfL?Gm9M@Hf@{+)-##@zZm7QdhSI{&OwDg;v5P8qbfi~v`_7oIsan6nZ|&2s
zllxb2+gz_?RQ)e^Js*^?L?&7OT*1J=z>K%tlw)9EKqM?kvkR2BjX~Z8x4(M(&vG3y
z5NLUC=Xzte?`j1l-)g~y<`cQ%oOk6mpXEOp_T??_Je`UUJ2o5NH%R?!B3zx&aVwJ1
zbsDqt>X_wk6&~E0cddHG>d?xAo<W<A1pYn|bj*!Wd7^!R+|>@ZV`nP2yu7?SL1n2i
z&mHg5>CB&e=e@7^k>;>9Vbd~Ufm!PquI;n=cet!}y{O#WS-erd0uNr5IXQFo%yS2q
zx?MYEpul;d<(YEs<|8cgzI{A;;1w^(=kHxEJYUDaz#vYR&%qIm>GQkW*ab6hb#v_P
z4EZ2+DOWnD&LimJcJqo&4-P&)cdq7e$l2iYZyfD3vMod+`FeV`>@`wu&o009ZsAJv
ztF38*lU03Au}(8mE0~c~!hS?S{gP73wTYhd4_4_WS^e$sJF#1W^R&(Vh5yqW_9kuG
z_WZ%^eGJ$3%KSaN?eGE>*|jkyGk+{*dMdQ&s=4j_y}rF`gzJ`^aCE&Szw(Tbu)*8I
zR}cIG`~2aRXeaW04lYeFeeV0((ZvTI=e`(m-u^68{o8>VJ>^WR7O}MX{k&(P`{35x
zYt|XBu0Cv<yd=kSF=yw~9&aYMMH`;mtQ1`9)5`ziQ~HOCUZ>Kp*m#sqm-@L(XS+?E
z@q##m$lKf<x$%sz;~4)Z{!TRD>W@1sCi{2AgDk%%Z_oNYe^IRbI^~K76VJlG9-$Lw
z9{6mxPa<MJI9qQ@+Ez~3=QxT{aFnM`IG6=)WmR`Q*u7XKck(tB%{<?Q*Eb$qcQ8=O
zIhXaergZ9iyY{))Wa~L=|F^$7ynLr>c}d*l+=zK;I*zMXT>o}!!}smR_cuBJ4*Fhq
zzkK26MwvLKSt|2xT72e}cAoK#A+9&WYW}4bGd=dAW!#3lv@aNBv`)H`bAe4`(;602
z_h~;3%%{0*>pRh%bHkaV@{IqYM@j0e@#)89#UlfY3iVE9ZC?B7YE<3kvX4J?YmOfL
z-^-m}_G|u(i>D?oS6tG?^J8PN#l$|l@1h!;6j!yb=S(#G&T~s@p4Wq)oNxI-8RNLu
zi;edg7#QsEW(>@-4_p|$jgBs!?IBioKi+ZU`5D%`u4F8;DcUo8mW5}p?ey&{m+jto
zgTeR3v1!>^|9;E8<Cws5V)mzNc|nCrzvdi%`ToxCsqgD<v^Rcn5q`HhL-}?1ES-xS
z#*f_lrGHLuk51M6|KacbeTS2JBULg~<l>U5`pZpRgf_U<^d5D)$0@izbn5f-KSN(U
zZYg}$seQ(2YURBCm020H7Po~cG_eFlUwRT0(D96;OIYl(%a19~&#!Css4*9va#~2^
z?{T+jO9ighZ;$&e*jl1|-(XJLhb6mz{Yz8MuJ&!eKWk-v{`T8zH{QA@wJdL@_t_g+
zr>Aa=eR2KV+pzba{zV<mES;XGalBK$yQXq>fXa%NB^`F&OG2KkI$hMC*u3=D(uGD(
zWll%k+ScrQ*sNoBgT40bjCt#VlcH~2OWAO1-elAK+gx&2B<$_&<85b$>r|c*xO(;e
zdpnV`so9S{ZuywHZ1oYBJ2y^TI}~|SLXs;^vM)(7cjbcCWDT)0pX(*LcMSp$C*-AE
z-Y7O>P0mrnmFpsw?Pg#*9a`e_;n2nneiL-PnN<5jQo^1a9{Lj0dQgfp`N3}`&&+=p
zKW92rq)B&ae3D5Lxa%@=+l`2|Vrp@g>t3d8o1}Up<Kg5Djeq2qZvA7B);VF3+v&iJ
z74A3P6N5KiZb;Lc@j3s=wY<}%3~GEY>^`45wfn_&kx#6$A^|(U$aY;w-d7MKz}3{(
z(^4gEdH7|Yko?`tPdIss!e=otd|&J@v7FOl<+a-f6yEDTdL;Dkb;ySKFHWaFb8k8?
zXfV_Az`a8sOKu;V5|LnV{ZbhB-7Plfk1Wfqc>H9>YxXVol^V;hI~usZ_rE>=`G&vE
zZYGYAtExXZTe%<Ftse19{n2h+<;^U2f~P#V`@HN%AJdz<b|J62&qv>h^J$h_-RAu<
zQ~u1SN&fHlYltZ@)NNX_M%8)J6}7x8bpb2aiigTaU$mO0?B92P%ajQd{7;<yE3=sE
zhF@%8l8xUC`AQwPr0JfqE9M6(9jM$M&?-D_=EiWFW1=gC&PGpbF~8!R$?)}2b{CU;
z7OUmF&-(ZOerzhRlsJ6#k;gQ*Y<s1W{S&7Ly#Ax$kmnF>I{AkU--Xl1FDsSVDVk($
zbQD_Z+0tZOb0L1Bf8d9pqOiO->2dd0*>XO;VDd;O*i-LO<WsJnzm$cybZ%x|=C@=H
z+v4D|EQhem3@<xpGI;fW_hq=ecDbq1E;HU)`xbEC_%ic8i^bssFP=BPG?KghJ^TI6
zB^($3G55+H^{aAn<YH;dHLuw6<GG{BBIT?LM(h`bGR`ZB?UVW17x;d~k~?X>_40SW
zD?WM5{^i7Pwo}uWRIPQq=ggSzEA)QhzJHn1cWKNt@><qAY31TAUw^)KSZJ+bRF|DU
z=lD#Oxwi}&ZiGZiaIcD(maeJoHNETgK^6w)rhOl}EX5{pykB{uV$0Q_?FNRTw<Cj&
zF4p!vIZwm;k;Ja&y0+z!JUcWFSx-yRnb)0jy4aUB{c)T7bgL;xTNfv!>^$eSecBBr
zOYV{s_KznI%FmB<(vn!_CsU~{`uHpJiu<|K9<cs?HTfRbInBElUS{Z8JvsLI&F8{j
zdHO;vtg4ro?{*s>j5_&0pY!tMWC!NS#=7VEJ!VaK!M;Pk$=Lk1cFM1J6D}XvJcG?Y
z{sK$MuQT~fa~c`f9ZXO1w$#XGa%ujPzbMke%II2zq_S1}L?L#8=I{U`uEWAb`zP-I
zTffHE^je>_*I#+5Z2tapv+g9%?O)!$NqKS+qmajhZ8{!@Gc0<T<}>%yTS@HRwNhWC
z;cjDJ{zJ{xQ(tQEUFz6sV)E+I(Rk*-YoSx#-zdH}C-2FtpT4I5Y=g^}tnA;G_;hk^
zp?{#`_ICb121cLjN_HnFH1PM-dOo(=<S#4xe^#mLgZcV}n#-Q#ZZFDwvPtBO<F2=H
z#|_-2_eI_HW51D;<rja$-a2Pb+L1Hy8@wg|H%x7B;a5M$?K7?O+*z|zA(CIKJJtyM
zmM5>Q6z4SOb(re^#3q1a(Q)JSHLWe#J&w#XvToT;&9*yM-FHfRoxXd9$@d>$rK@v_
zbMBV#Df1;q*=fsmNqoC&cJ{VRre8Vd-lECAn>RdibaZKY9?LLgseZtos_XYSLw|jr
zC%l{Qzs<%+?lA{q;_k1TyezrnvdW1_j<D0Z#i9EzU2QY7W_{4i^V7@s<bsg-n){{-
zuQ;RnJU{Wda^)_~i2+}aHZV7tWKG<`csTym(*4RSmOTH}D;ywZ{p)~%R^D%KzDc_p
z<x_7xuX(yJ^oDWx{w9$dOIK!1TP3ORnJcbt-%ClO^VxfBr<X1`lz8vNkuAp;D)HyZ
z|9hdc;_sSxuNRl+?c-Cc=-sUEyWz<j!Sf9H`{p}2?9Y@?PI&8lTFoeb=C^-b>#l{z
z`4#YAoopqW{r63@TD_)tjZg8z?gwYGXHSn{KAUr{-p9)S^Ugf6WlMf~o#6@CFLJzX
zdH<hZ9df~H4i}TTPNwx$9$YOK`Kf53TI+^S-FEAhv={pS<4FuvaT63?v1azo<}3^C
z`JV1|wr^_ULXtT5oH~Eu;iMamPnIwA7VfMRa`*qxqZ+hSaX*up)<O%uKT3t33M^jb
zJlS?83PM*{HG3B_G_)zKtz-S9utT*g;oPd0Sl(4(A>Yo)om&0no8S5SE-B?%Jx`WC
zKJMLHH(|k{pjUrHiurGQJa0MDv0h}$tVM?>PY$=5G^1#f(QED4^n3M%(t!@e8Pm@T
zoD*m`8dUzna!zfPxfJ{1^9NR3PR=-MC47SYRM1PoxtdGfi?A0OrIlFoG#4+vW%289
z;2V|W{pPo`&q;pUZ1+G`q~xQ~c@u>dDe}g@`3`W-&uU+_R7@#*{km@jT^}sY{^b@o
zsfp)fUY@BPpul{sMVaRwht`9)uj<ce9RKOIndeQOPImd`#j($Y+BGZG-AfLJ<$g|o
zb+P5y;=3$qi~N<%mra~~`NZ1HC*SH7jCpS)P2rX?JtmU3$bX5s+NIgTFV}W{`R4KK
zkzHNbhLkmn`NiT7E58yy7gk>5@%+8istfBkc5>c}d1>gf=9+HYx5O9w8*2W8yGK`4
z7rZ^q!NA~Pgl||0M}raEAlf=5I{&tU$npE%HI?oeZtT@xIV;3d)_ZDL<}It0>m~Yx
zO_;?Vq+Q<f-tNt_)V)hXcf0e=x&L#%^mB3hQ1ji)ua{i(7I7~&VLsbCp=!I!{rs$$
zMcr4I7e745^7%|<^21Zjvsle#bMiXB-TR*2VzH0?XiQ4gg@>C=!wlv7%G`Z!d+|+u
z*v$0S(`f$l7XHSid7<rTjk?UARz)qk9~yeKWbewtdDk?TMMph4Yj*mBn%?H3*!8D#
zgZ?|N4m@9bVEu(%_5I5E8$L}mP~B;IZj0R`H(ArW*?|vk&X9fR_ItsHXW{%EH4<C5
zFI2ncT`)QK?cS+p9~_+9RcIlz$k}SDfk1Lm_1~YjzkJ)@Q7z`xT{E$8TEae$BA=Uk
zj;uK$>{B_<EOYWv&tH2MnJi}&@76I9N%Fb8jN{yvgBu0*oY3rNJ$-L}pT)Urrk8Wr
zFHiY-;mBmsx`{fa*;UWJCwDZ}BwowjFg?*_&C597x+K}EWA!Xk`_s-d)w!DY{i)k<
z*F|`agZa`o9x_Lt==RC2DN?tcr#Fdbb&~1pj?5z-3jAS{Z2lO9L52?{f^#=6Vq{>*
z!8=$EYV$@r=jWBA=9R>UR2HNb$Ldw&<{S+>op;DUV9#gKz;&XA{3};!F}2(E=?1O5
z68%)o`;m}+iZoaK{jHakx*8lNU-Dj4QJz+2maW}<agz8!rY(s<M-t6nwJ>(LEb~bz
z{9&^*KfZNwVzxn~rUBzKF5&i_TVKr*nEd3x%FZQBjT2X1nG;oWRWi_I=0zTl^_O?5
zbc%?z<VGJC<YV7%u=MjGwsP+~dZmJLZ04!+CU-fn3OOj0W#|(#<;{^yjWdi7oR&mi
z-MU|^Ti|7}&$%f(Z>(S7vQ72Cj^bzE%9ZblzG?fxkR-S9?yRK(JuA2IPf2Eq6nk=Y
zrlNe%`Fl&|xJCqB*|694(CkUeSDt)1MRmfdlyibYl1m$xsj+W}*b}epH*>{>?G|nU
zcDuJ`e=rekI<}6>!1wmoCGJ~pG4p)%uHA3n9q~E5RAIT)CH2)GXPDpX@O|f8I)%Bu
zPw}|d@hGk4@QH1g-`2ghEQwYW%s9<g|9Y!Z(DirvJ5H;%ehGP-lDFo`)~DXHl#iZ%
zA|5A@));>F`$d;~0e>%VQj`6E$Y(vL@zV>|yZ-H5_apYv-v1)+=f&69I)}=C>`z~D
z|3%}9q}9pGOcnku&;P**>gGJZc4p2dMg|6DE>PBHWD;ROtgA*|Neo(7jkL@fGK~uv
z0s*boL0M%T;Eiey^2$=sx?_Yn8lah(9F%q1=$eri)qzqTLURH$c&Rsd5h~JBRCEK7
zR~~`Z-XRQF#0EA1u@VnmEApH+sB}YU-Ny;m3i3L5g%a{YC3I7eXYxSv-w0FY^Mg&n
zm<>la0eKP&G~<Ub;k^*b{2#hT<ard(92Y{PnFy8{FLdq5(+Qx-7KHXUqDb1IQ!nVc
zk%!7rV@^dJNjGFf4OD`l8-P6Ai)uid1oB86XbcWK9)oTG^7tic{4q<R8-N%iMmGg{
za1hlQ#nR}efJaBr4S0_<#s(TgL<D1u46*@`>;&=+xXuKRE}@%+JZgn%R+ubyv%pTm
zY8LW15UN=La@frRhcZ^PkcUf9gW5|T*({Lnz(Eam7P{|{8ycuVov(my3OJ0=4M48J
sP@Pe&iEaSc8CXq0E)-Fnk*v+YfTivW@MdKL#VV-uv|?po2-5}e0JZ2MSO5S3

literal 0
HcmV?d00001

diff --git a/scenarioProcessing/relativizeShocks.R b/scenarioProcessing/relativizeShocks.R
index 6a2fdc1..e012ebd 100644
--- a/scenarioProcessing/relativizeShocks.R
+++ b/scenarioProcessing/relativizeShocks.R
@@ -2,11 +2,15 @@
 # Uses the stocks in the flood scenario 2 and in the country level aggregate stocks derived
 # from it to relativize the shocks in the earthquake and flood1 scenarios.
 # 
+# Benjamin Blanz 2024
+# 
 
 source('funRelData.R')
 
-nuts3LevelStocks <- read.csv("helperData/nuts3LevelStocks.csv", row.names=NULL)
-countryLevelStocks <- read.csv("helperData/countryLevelStocks.csv", row.names=NULL)
+nuts3LevelStocksNACE <- read.csv("helperData/nuts3LevelStocksNACE.csv", row.names=NULL)
+nuts3LevelStocksGTAP <- read.csv("helperData/nuts3LevelStocksGTAP.csv", row.names=NULL)
+countryLevelStocksNACE <- read.csv("helperData/countryLevelStocksNACE.csv", row.names=NULL)
+countryLevelStocksGTAP <- read.csv("helperData/countryLevelStocksGTAP.csv", row.names=NULL)
 
 
 files <- list.files('scenarios',pattern = 'csv',recursive = T)
@@ -16,8 +20,13 @@ for(f.i in 1:length(files)){
 	file <- files[f.i]
 	cat(sprintf('%i of %i %s\n',f.i, length(files),file))
 	data <- read.csv(file,row.names=NULL)
-	data.CNT <- relData(data,nuts3LevelStocks,countryLevelStocks)
-	write.csv(data.CNT,gsub('.csv','-rel.csv',file),row.names = F)
+	if(grepl('GTAP',file)){
+		data.rel <- relData(data,nuts3LevelStocksGTAP,countryLevelStocksGTAP)
+		write.csv(data.rel,gsub('.csv','-rel.csv',file),row.names = F)
+	}else{
+		data.rel <- relData(data,nuts3LevelStocksNACE,countryLevelStocksNACE)
+		write.csv(data.rel,gsub('.csv','-rel.csv',file),row.names = F)
+	}
 }
 
 
diff --git a/scenarioProcessing/stockAggregation.R b/scenarioProcessing/stockAggregation.R
index 9c7cdad..37e4a4e 100644
--- a/scenarioProcessing/stockAggregation.R
+++ b/scenarioProcessing/stockAggregation.R
@@ -1,8 +1,13 @@
 # 
-# Aggregate NUTS3 scenario data to country level
+# Aggregate NUTS3 scenario stock data to country level
 # 
 # library(sf)
-# worldAdmin1 <- read_sf('scenarios/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp')
+# worldAdmin1 <- read_sf('scenarios/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp') 
+# 
+# Benjamin Blanz 2024
+# 
+
+sectorColPattern <- 'ALL|TOTAL|^[A-Z]$|AGR|MIN|MFG|EGW|CNS|TRD|OTP|WTP|CMN|OFI|OBS|REA|PUB|OSG'
 
 # read scenario file with NUTS3 data ####
 library(readxl)
@@ -37,7 +42,7 @@ sink('helperData/nuts3LevelStocksMetadata.csv')
 cat(sprintf('Unit %s,,\n',stockUnit))
 cat(sprintf('Labels,, \n'))
 cat(sprintf('Sector,Label,Type\n'))
-sectorCols <- grep('ALL|TOTAL|^[A-Z]$',names(stocksNUTS3),perl = T)
+sectorCols <- grep(sectorColPattern,names(stocksNUTS3),perl = T)
 for(n in colnames(stocksNUTS3)[sectorCols]){	
 	cat(sprintf('"%s", "%s", "%s"\n',n,stockLabels[n],stockTypes[n]))
 }
@@ -45,7 +50,7 @@ sink()
 
 
 # aggregate to country level
-sectorCols <- grep('ALL|TOTAL|^[A-Z]$',names(stocksNUTS3),perl = T)
+sectorCols <- grep(sectorColPattern,names(stocksNUTS3),perl = T)
 stocksCNT <- aggregate(stocksNUTS3[,sectorCols],
 											 by = list(Category=stocksNUTS3$CNTR_CODE),FUN=sum)
 colnames(stocksCNT)[1] <- 'CNTR_CODE'
@@ -65,7 +70,7 @@ sink('helperData/countryLevelStocksMetadata.csv')
 cat(sprintf('Unit %s,,\n',stockUnit))
 cat(sprintf('Labels,, \n'))
 cat(sprintf('Sector,Label,Type\n'))
-sectorCols <- grep('ALL|TOTAL|^[A-Z]$',names(stocksNUTS3),perl = T)
+sectorCols <- grep(sectorColPattern,names(stocksNUTS3),perl = T)
 for(n in colnames(stocksNUTS3)[sectorCols]){	
 	cat(sprintf('"%s", "%s", "%s"\n',n,stockLabels[n],stockTypes[n]))
 }
-- 
GitLab