Add documentation and GTAP capability

.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.*#

scenarioProcessing/a1ProcessScenarioData.R

+#
+# 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')

scenarioProcessing/aggregateShocks.R

 #
 # aggregate the flood2010 and earthquake scenarios to country level
 #
+# Benjamin Blanz 2024
+#
 
 source('funAggregateNuts2CNT.R')
 codes <- read.csv("helperData/nuts3fid4Codes.csv")

scenarioProcessing/funAggregateNuts2CNT.R

+#
+# 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'

scenarioProcessing/funRelData.R

+# 
+# 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])){

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

scenarioProcessing/stockAggregation.R

 # 
-# 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') 
+# 
+# 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]))
 }