diff --git a/CommonFunctions.R b/CommonFunctions.R
index 6283d01d422b93cbf4992befee5e9aa48d470faa..737d3293ea8fa02fd6a4918b346416277a02c98a 100644
--- a/CommonFunctions.R
+++ b/CommonFunctions.R
@@ -4,7 +4,6 @@ require(reshape2)
 data_dir_root = "~/Documents/LURG/Data"
 fao_data_dir = file.path(data_dir_root, "FAOStat-April2019")
 
-
 getItemDetails = function() {
 	itemDetails = fread("~/Documents/R_Workspace/landalloc/itemsInputs.csv")
 	itemDetails[, allowedForVeggie:=TRUE]
diff --git a/PLUM/Cluster/MonteCarloStuff.R b/PLUM/Cluster/MonteCarloStuff.R
index 619147bfd1d84683057a112e4b637cf81d8e34ec..e8092841b8bb4c9f99e2f9721d2b8ba80ba70e07 100644
--- a/PLUM/Cluster/MonteCarloStuff.R
+++ b/PLUM/Cluster/MonteCarloStuff.R
@@ -83,7 +83,7 @@ generateMCSample <- function(n, vals, ensem) {
 #paramDt = rbind(paramDt, data.table(scenario="hind", pname="NORMAL_WORKS_TOO", dist="norm", mn=5, sd=2), fill=TRUE)
 #write.table(paramDt, file.path(baseSimDir, "hindparms.csv"),sep=",",row.names=FALSE,quote=FALSE)
 
-baseSimDir="/Users/peteralexander/Documents/R_Workspace/UNPLUM/data/sims"
+baseSimDir="C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2\\data\\sims"
 
 generateScenarios = function(inputFile, outputFile, numScen, drawPlot=FALSE, ensemblePrefix="") {
 	paramDt = fread(file.path(baseSimDir, inputFile))
@@ -140,9 +140,13 @@ generateScenarios = function(inputFile, outputFile, numScen, drawPlot=FALSE, ens
 }
 
 generateScenarios("ssp_params.csv", "ssp_sims11.csv", 50, ensemblePrefix="ssp11/")
-generateScenarios("rcp_params.csv", "rcp_sims.csv", 50)
+generateScenarios("rcp_params_climate.csv", "rcp_sims_climate_ns.csv", 5, ensemblePrefix="climateShocks/")
 generateScenarios("hind_params_ha.csv", "hind_sims_ha.csv", 50)
 
-generateScenarios("shocks_params.csv", "shocks_sims.csv", 10, ensemblePrefix="shocks/")
+generateScenarios("shocks_params_ssp2.csv", "baseline.csv", 5, ensemblePrefix="shocks/baseline/SSP2")
 
+generateScenarios("shocks_params_ssp2.csv", "demandSystem.csv", 5, ensemblePrefix="demandSystem/baseline/")
 
+generateScenarios("prot_params.csv", "protectionismLevels.csv", 5, ensemblePrefix="protectionismLevels/")
+
+generateScenarios("prot_params.csv", "protectionismLevels.csv", 5, ensemblePrefix="protectionismLevels/")
diff --git a/PLUM/SetupData/CreateYieldShockMaps.R b/PLUM/SetupData/CreateYieldShockMaps.R
index fedf9e89327045505c669ccce35d421487def5c3..b111d56758eb9f7a5bd50f4788eba2bce38a766a 100644
--- a/PLUM/SetupData/CreateYieldShockMaps.R
+++ b/PLUM/SetupData/CreateYieldShockMaps.R
@@ -1,15 +1,52 @@
 require(data.table)
 require(raster)
 
-data_dir_root = "~/Documents/R_Workspace/UNPLUM/data"
+#data_dir_root = "~/Documents/R_Workspace/UNPLUM/data"
+data_dir_root = "C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2\\data"
 
 createMapForCountries = function(countries, outputFileName) {
 	country_codes = fread(file.path(data_dir_root, "country_codes4.csv"), header=TRUE)
 	boundary = raster(file.path(data_dir_root, "halfdeg", "country_boundaries.asc"))
 	crs(boundary) <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
 	
-	b2 = boundary %in% country_codes[Country %in% countries, numCode]
+#	b2 = boundary %in% country_codes[Country %in% countries, numCode]
+	b2 = boundary %in% country_codes[Region %in% countries, numCode]
 	writeRaster(b2, filename=file.path(data_dir_root, "halfdeg", "yieldshockmaps", outputFileName), format='ascii', bylayer=TRUE, overwrite=TRUE)
 }
 
-createMapForCountries(c("China", "Brazil"), "china_and_brazil.asc")
\ No newline at end of file
+ga = createMapForCountries(c("Latin America & Caribbean"), "Latin America & Caribbean.asc")
+
+
+#year,mapFilename,crop,shockFactor
+climateShockProbs = fread(file.path("C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2\\data","climateShockProbabilities.csv"))
+scenarioTable = fread(file.path("C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2\\data\\sims\\cyberauto.csv"))
+climateMaps = basename(list.files(path = "C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2\\data\\halfdeg\\yieldshockmaps", full.names = TRUE, recursive = TRUE))
+
+
+
+shocksHappen = data.table(SSP_Scenario = scenario[,SSP_SCENARIO], Scenario = scenario[,Scenario], Year = seq(2010,2100,1), 
+		climate =as.logical(rbinom(91,size=1,prob = climateShockProbs[rcp == scenario[,RCP],climate])))
+
+tab = NULL
+for(year in shocksHappen[,Year]){
+		if(shocksHappen[Year == year, climate]){
+
+			mapFilename = sample(climateMaps, 1)
+			shockFactor = runif(1,0.1,0.5)
+			tab = rbind(tab, cbind(year,mapFilename,crop='maize',shockFactor))
+			tab = rbind(tab, cbind(year,mapFilename,crop='rice',shockFactor))
+			tab = rbind(tab, cbind(year,mapFilename,crop='cereal',shockFactor))
+			}		
+		
+}
+
+
+for(i in 1:nrow(scenarioTable)){
+#first if column exists
+	if("SHOCKS_POSSIBLE" %in% colnames(scenarioTable)){
+		if(scenarioTable[i,SHOCKS_POSSIBLE]){
+			print(paste("Shocks possible in ",scenarioTable[i,Ensemble],scenarioTable[i,Scenario]))
+			createShockFile(scenario = scenarioTable[i])
+		}
+	}
+}
\ No newline at end of file
diff --git a/PLUM/SetupData/NutritionDataExtract.R b/PLUM/SetupData/NutritionDataExtract.R
index 18e9e390c2c7f59d50f0e07ed5a2c4711de8f096..c6ae08ab3203290395048990005b9a9bc6d8b787 100644
--- a/PLUM/SetupData/NutritionDataExtract.R
+++ b/PLUM/SetupData/NutritionDataExtract.R
@@ -3,7 +3,7 @@ processedOilCrops = c("Rape and Mustard Oil", "Rape and Mustard Cake", "Soyabean
 			"Soyabean Cake", "Sunflowerseed Oil", "Sunflowerseed Cake", "Palmkernel Oil")
 
 getItemGroupMapping = function () {
-	faoItemGroupMapping = fread(file.path(fao_data_dir, "ItemGroup.csv"))
+	faoItemGroupMapping = fread(file.path(fao_data_dir, "nutritionItemGroups.csv"))
 	setnames(faoItemGroupMapping, gsub(' ', '', (names(faoItemGroupMapping))))  
 	faoItemGroupMapping[, Item := gsub(',', '', (Item))]
 	
@@ -14,7 +14,7 @@ getItemGroupMapping = function () {
 					"Wheat and products", "Barley and products", "Oats", "Maize and products", "Millet and products", "Sorghum and products", "Rice (Paddy Equivalent)"),	
 			fcr=c(2.3, 3.3, 6.4, 0.7, 25, 15, 15, 9.5, rep(NA, 7)))
 	
-	itemGroupMapping = rbind(itemGroupMapping, faoItemGroupMapping[ItemGroup %in% c("Oilcrops", "Pulses", "Starchy Roots"), list(plumCropItem=ItemGroup, Item, fcr=NA)])
+	itemGroupMapping = rbind(itemGroupMapping, faoItemGroupMapping[plumCropItem %in% c("Oilcrops", "Pulses", "Starchy Roots"), list(plumCropItem, Item, fcr=NA)])
 
 	itemGroupMapping = rbind(itemGroupMapping, data.table(plumCropItem = rep("Oilcrops",length(processedOilCrops)), Item=processedOilCrops, fcr=NA))
 	setkey(itemGroupMapping, plumCropItem)
diff --git a/PLUM/SetupData/PreparePLUMData.R b/PLUM/SetupData/PreparePLUMData.R
index 56121784e84d0b33f0199d14c884c3d7b96aa6f7..4652c3953ca7157631cb0dbbc83819cf4bf5528b 100644
--- a/PLUM/SetupData/PreparePLUMData.R
+++ b/PLUM/SetupData/PreparePLUMData.R
@@ -1,6 +1,6 @@
 require(raster)
 #write files for java model
-#output_dir_root = "C:\\Users\\rhenry2\\Desktop\\PLUM_input"
+#output_dir_root = "C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2"
 output_dir_root = "/Users/peteralexander/Documents/R_Workspace/UNPLUM"
 
 # Stock levels
@@ -20,71 +20,74 @@ writePlumFiles = function(dt, year, output_dir) {
 				gdpRatio=c(15.665/3.874, NA, 58.641352/115.419051, 9.716103/3.140509))
 	} else{
 		# From country_data in 2010, supplemented by http://data.worldbank.org/indicator/SP.POP.TOTL?locations=SD
-		missingCountries = data.table(Country=c("Chad", "Chad", "Chad", "Thailand", "Malawi", "Algeria"), 
-				toCountry=c('Sudan', 'South Sudan', 'Somalia', 'Myanmar', "Democratic Republic of the Congo", "Libya"), 
-				popRatio=c(36.115/12.29851, 10.056/12.29851, 9.581/12.29851, 51.73301/66.69202, 65.93871/14.76982, 6.278438/36.03616), 
-				gdpRatio=c(65.634/10.658, 15.727/10.658, 5.352/10.658, 49.541/210.090, 15.665/3.874, NA)) # 2013 for Somalia GDP
+		missingCountries = data.table(Country=c("Chad", "Chad", "Chad", "Malawi", "Algeria"), 
+				toCountry=c('Sudan', 'South Sudan', 'Somalia', "Democratic Republic of the Congo", "Libya"), 
+				popRatio=c(36.115/12.29851, 10.056/12.29851, 9.581/12.29851, 65.93871/14.76982, 6.278438/36.03616), 
+				gdpRatio=c(65.634/10.658, 15.727/10.658, 5.352/10.658, 15.665/3.874, NA)) # 2013 for Somalia GDP
 	}
 	
 	#SSP
 	newSsp = rbind(dt$ssp[!Country %in% missingCountries[!is.na(gdpRatio)]$toCountry], 
 			merge(missingCountries, dt$ssp, by="Country")[, list(Model, Scenario, Year, population=population*popRatio, Country=toCountry, gdp=gdp*gdpRatio)])
 	write.table(newSsp[!is.na(population) & !is.na(gdp) & Year >= year, list(Model,Scenario,Year,population,Country,gdp)], file.path(output_dir,"ssp.csv"), sep=",", row.names=FALSE, quote=FALSE)
-	
+
 	# Base Consumption (of commoditites)
 	newBaseConsumption = dt$baseConsumption[!is.na(baseCpc), list(Country, basePop, Item, baseCpc)]
 	newBaseConsumption = rbind(newBaseConsumption, merge(missingCountries, newBaseConsumption, by="Country")[, list(Country=toCountry, basePop=basePop*popRatio, Item, baseCpc=baseCpc)])
-	popForG = processedFull$histDt[Year == 2010 & Item == "Cereals" & !is.na(population), list(Country, popForGrouping=population)]
+	popForG = dt$histDt[Year == 2010 & Item == "Cereals" & !is.na(population), list(Country, popForGrouping=population)]
 	newBaseConsumption = merge(newBaseConsumption, popForG, by='Country', all.x=TRUE)
 	newBaseConsumption[is.na(popForGrouping), popForGrouping:=basePop]
 	write.table(newBaseConsumption[!Country %in% c("World", "Sudan (former)")], file.path(output_dir,"base_consump.csv"), sep=",", row.names=FALSE, quote=FALSE)
-	
+
 	# Net imports and production (of crops)
 	consumpDt = dt$cb_fs[Year == year]
 	consumpDt[!is.na(fcr), c('net_import', 'prod') := list(net_import * fcr, prod *fcr)]
 	consumpAggreated = consumpDt[, list(net_import=sum(net_import, na.rm=TRUE), prod=sum(prod, na.rm=TRUE)), by=list(plumCropItem, Country, Year)]
 	consumpAggreated = rbind(consumpAggreated, merge(missingCountries[!toCountry %in% consumpDt$Country], consumpAggreated, by="Country")[, list(Country=toCountry, Year, plumCropItem, net_import=net_import*popRatio, prod=prod*popRatio)])
 	write.table(consumpAggreated, file.path(output_dir, "net_imports.csv"),sep=",",row.names=FALSE,quote=FALSE)
-	
+
 	# Animal number mapping
 	animalNum = merge(consumpAggreated, dt$livestockNumbers[Year == year, list(Country, plumCropItem=plumDemandItem, LivestockItem, animalM)], by=c("Country", "plumCropItem"))
 	write.table(animalNum[, list(Country, plumCropItem, LivestockItem, animalRate=animalM/prod, prod)], file.path(output_dir, "animal_numbers.csv"),sep=",",row.names=FALSE,quote=FALSE)
-	
+
 	#stock levels
 	cb = copy(dt$cb_fs)
 	cb[, StockLevel := calculateStockLevels(stockvar), by=c("Country","Item")]
 	globalstock = cb[Year == year, list(StockLevel=sum(StockLevel)), by=c("plumCropItem", "Item", "fcr")]
 	st = rbind(globalstock[is.na(fcr), list(sum(StockLevel)), by=plumCropItem], globalstock[!is.na(fcr), list(sum(StockLevel*fcr)), by=plumCropItem])
 	write.table(st, file.path(output_dir, "global_stocks.csv"), sep=",", row.names=FALSE, quote=FALSE)
-	
+
 	# Bioenergy demand
 	bio = dt$cb_fs[Year == year & !Item %in% itemDetails[animal==TRUE], list(Country, Year, Item, other=other, Item, plumDemandItem)]
 	
 	if (exists("dt$isHindcast") && dt$isHindcast){
-		dt = bio[, list(other=0), by=list(Country, Year, Item=plumDemandItem)] ##set to zero for hindcast
+		bio = bio[, list(other=0), by=list(Country, Year, Item=plumDemandItem)] ##set to zero for hindcast
 	} else {
-		dt = bio[, list(other=sum(other, na.rm=TRUE)), by=list(Country, Year, Item=plumDemandItem)]
+		bio = bio[, list(other=sum(other, na.rm=TRUE)), by=list(Country, Year, Item=plumDemandItem)]
 	}
-	write.table(dt, file.path(output_dir, "bio_demand.csv"), sep=",", row.names=FALSE, quote=FALSE)
+	write.table(bio, file.path(output_dir, "bio_demand.csv"), sep=",", row.names=FALSE, quote=FALSE)
+
+	#energy per mass values 
+	cals = merge(dt$cb_fs[Year == year],dt$country_data[Year == year, list(Country,population, gdp_pc)], by=c('Country'))
+	cals[!is.na(fcr), c('food') := list(food * fcr)]
+	calsB = cals[,list(gdp_pc=mean(gdp_pc),kcalPerT=sum(energy_kcal_pc*population*365,na.rm=TRUE)/sum(food)),by=c('Country','plumDemandItem')]
+	calsB = rbind(calsB[,list(Country,plumDemandItem,kcalPerT)], merge(missingCountries, calsB, by="Country")[, list(Country=toCountry, plumDemandItem, kcalPerT)])
+	calsC = calsB[!Country %in% c("World", "Sudan (former)") & kcalPerT != 0]
+	write.table(calsC, file.path(output_dir,"calories_per_t.csv"), sep=",", row.names=FALSE, quote=FALSE)
 	
 	# base_demand_fracts
-	newComBal = rbind(processedFull$cb_fs[,list(Country,Year,Item,supply=(prod+net_import-stockvar),net_import,food,processing,feed,other)], 
-			merge(missingCountries, processedFull$cb_fs, by="Country")[, list(Country=toCountry, Item, Year,
-							supply=(prod+net_import-stockvar)*popRatio, net_import=net_import*popRatio,
-							food=food*popRatio, processing=processing*popRatio, feed=feed*popRatio, other=other*popRatio)])
+	newComBal = rbind(dt$cb_fs[,list(Country,Year,Item,net_import,food,feed,other)], 
+			merge(missingCountries, dt$cb_fs, by="Country")[, list(Country=toCountry, Item, Year,
+							net_import=net_import*popRatio, food=food*popRatio,feed=feed*popRatio, other=other*popRatio)])
 	
 	cereals = merge ( itemDetails[plumDemandItem %in% c("Cereals", "OilcropsPulses"), list(Item, plumDemandItem, plumCropItem)],
-			newComBal[Year == 2011, 
-					list(Country, Year, Item,
-							supply=supply, 
-							net_import=net_import, 
-							food_rate=(food+processing)/supply,
-							feed_rate=feed/supply,
-							other_rate=other/supply)],
-			by = "Item")
+			newComBal[Year == 2010,	list(Country, Year, Item,net_import,food,feed,other)],by = "Item")
 	
-	cereals2 = cereals[, list(food=sum(supply*food_rate, na.rm=TRUE)), by=list(Country, plumDemandItem, Item=plumCropItem)]
-	write.table(cereals2, "~/Documents/R_Workspace/UNPLUM/data/base_demand_fracts.csv", sep=",", row.names=FALSE, quote=FALSE)
+	cereals2 = cereals[, list(food=sum(food, na.rm=TRUE)), by=list(Country, plumDemandItem, Item=plumCropItem)]
+	write.table(cereals2, file.path(output_dir,"base_demand_fracts.csv"), sep=",", row.names=FALSE, quote=FALSE)
+
+
+
 }
 
 writePlumFiles(dt=processedFull, year=2010, output_dir=file.path(output_dir_root, "data"))
@@ -139,7 +142,7 @@ missing_owu = missing_owu[, list(FPU, Year, Value=value2001*(1+rate)^(Year-2001)
 write.table(rbind(missing_owu, owu[Year>2000]), "~/Documents/R_Workspace/UNPLUM/hind1970/data/other_water_uses.csv", sep=",", row.names=FALSE, quote=FALSE)
 write.table(rbind(missing_owu, owu[Year>2000]), "~/Documents/R_Workspace/UNPLUM/data/other_water_uses.csv", sep=",", row.names=FALSE, quote=FALSE)
 
-# FEED CHECKS/CALCS
+# FEED CHECKS/CALCS TODO feeds are in raw equiv terms for processed goods so DM values need altered to match 
 fbs = fread("/Users/peteralexander/Documents/R_Workspace/UNPLUM/output/fbs.txt")
 fao = processedFull$cb_fs[Country=="World" & Year == 2010, list(faoFeed=sum(feed), faoFeedDM=sum(feed*dryMatter)), by=list(Crop=plumCropItem)]
 #plum=merge(fbs, itemGroupMapping[, list(drymatter=mean(drymatter)), by=list(Crop=gamsCrop)], by="Crop")[Year == 2010, list(mono=sum(MonogastricsFeed), rum=sum(RuminantsFeed), monoDM=sum(MonogastricsFeed * drymatter), rumDM=sum(RuminantsFeed * drymatter)), by=Crop]
@@ -158,12 +161,6 @@ fao = processedFull$cb_fs[Country=="World" & Year == 2010, list(faoFeed=sum(feed
 # It also implies we are missing 495 Mt of feed, primarily from forage crops
 # This is 12.7% of total feed requirements, 495 / (1162 + 2754)
 
-#energy per mass values 
-
-cals = merge(processedFull$cb_fs[Year == 2010],processedFull$country_data[Year == 2010, list(Country,population, gdp_pc)], by=c('Country'))
-cals[!is.na(fcr), c('food') := list(food * fcr)]
 
-calsB = cals[,list(gdp_pc=mean(gdp_pc),kcalPerG=sum(energy_kcal_pc*population*365/1000000,na.rm=TRUE)/sum(food)),by=c('Country','plumDemandItem')]
 
-calsB[Country=='World']
 
diff --git a/PLUM/SetupData/ReadFaoForPlum-notHindcast.R b/PLUM/SetupData/ReadFaoForPlum-notHindcast.R
index cf05403e010b460d6d3d6bcb20bc63871f4c88f9..e1ccdd43680987643ae87541af31ba29daf03fb2 100644
--- a/PLUM/SetupData/ReadFaoForPlum-notHindcast.R
+++ b/PLUM/SetupData/ReadFaoForPlum-notHindcast.R
@@ -79,40 +79,37 @@ disaggregateCountries = function(Dt,wbCountryDt, colsToChange){
 	Dt
 }
 
+
 getHistoricDataset  = function(commodityDt, countryDt) {
 	
 	h_dt = commodityDt[, list(Country, Year, Item, plumDemandItem, dryMatter, fcr,
-					supply=prod+net_import-stockvar, 
 					net_import=net_import, 
-					food_rate=(food)/(prod+net_import-stockvar),
-					feed_rate=feed/(prod+net_import-stockvar),
-					other_rate=other/(prod+net_import-stockvar))]
+					food=food,
+					feed=feed,
+					other=other)]
 	
-	h_dt[, dm_food := supply*food_rate*dryMatter]
+	h_dt[, dm_food := food*dryMatter] ##TODO this doesn't work now because processed goods have been converted to their raw equivelents
 	
 	# animal product feed conversion adjustment
-	h_dt[!is.na(fcr), c('supply', 'net_import') := list(supply * fcr, net_import * fcr)]
+	h_dt[!is.na(fcr), c('food', 'net_import','feed','other') := list(food * fcr, net_import * fcr, feed * fcr, other *fcr)]
 	
-	h_dt2 = h_dt[, list(supply=sum(supply), net_import=sum(net_import), 
-					food_rate=sum(supply*food_rate, na.rm=TRUE)/sum(supply, na.rm=TRUE),
-					feed_rate=sum(supply*feed_rate, na.rm=TRUE)/sum(supply, na.rm=TRUE),
-					other_rate=sum(supply*other_rate, na.rm=TRUE)/sum(supply, na.rm=TRUE),
-					avg_fcr_dm_adj=sum(dm_food, na.rm=TRUE)/sum(supply*food_rate, na.rm=TRUE) # factor to convert from wet to dry, and also from feed equivalent quantities back to real quantities
+	h_dt2 = h_dt[, list(food=sum(food, na.rm=TRUE),
+					net_import=sum(net_import, na.rm=TRUE), 
+					feed=sum(feed, na.rm=TRUE),
+					other=sum(other, na.rm=TRUE),
+					avg_fcr_dm_adj=sum(dm_food, na.rm=TRUE)/sum(food, na.rm=TRUE) #TODO factor to convert from wet to dry, and also from feed equivalent quantities back to real quantities
 			), by=list(Item=plumDemandItem, Country, Year)]
 	
 	h_dt2 = merge(h_dt2, countryDt[, list(Country, Year, gdp,population)], by=c('Country', 'Year'))
-	h_dt2[, c('gdpPc', 'cpc') := list(gdp/population, food_rate*supply/population)]
+	h_dt2[, c('gdpPc', 'cpc') := list(gdp/population, food/population)]
 }
 
 getBaseConsumption  = function(year,consumpDt, gdpDt) {
 	dt = consumpDt[Year == year, list(Item, Country, baseYear=Year, 
-					baseSupply=supply,
 					avg_fcr_dm_adj,
-					baseFoodRate=food_rate,
+					baseCpc=cpc,
 					basePop=population, 
 					baseGdpPc=gdp/population)]
-	
-	dt[, baseCpc := baseSupply/basePop*baseFoodRate]
 	dt
 }
 
@@ -185,59 +182,52 @@ getProjectAllConsumption  = function(consumptionDt,sspDt, comCurvesDt) {
 	dt
 }
 
+extractionValues = function(dt, dtProcessed){
+
+	# Get commodities produced from processing
+	dt1= merge(dtProcessed[Year==2010,list(prod=sum(prod)),by=c('Country','processed')],
+			dt[Year==2010,list(Country,processed=Item,rawAmountProcessed = processing)],
+			by=c('Country','processed'))
+	
+	dt1 = dt1[prod > 0 & rawAmountProcessed > 0]
+	dt1[,extractionRate:=prod/rawAmountProcessed]
+
+	
+	extractionValues = dt1[Country=='World',list(processed,extractionRate)]
+	
+	extractionValues
+}
+	
 processData = function(minYear=1960, year=2010){
 	print(paste("Base year is ", year))
 	
-	itemMapping = fread("~/Documents/R_Workspace/landalloc/itemsProcessingMap.csv")
 	com_bal = readCommoditiyBalance()
 	food_supply = readFoodSupply()
-	
+	cropProduction = readCropProd()
 	# Original FAO data for commodities of interest
 	cb_fs = merge(com_bal, food_supply, by=c("Country", "Year", "Item"), all.x=TRUE)
 	cb_fs = merge(cb_fs, itemDetails[, list(Item, plumCropItem, plumDemandItem, fcr, dryMatter)], by="Item")
 	
+	
 	cb_fs_sugar = cb_fs[Item %in% c("Sugar cane", "Sugar beet"), lapply(.SD, sum, na.rm=TRUE), by=list(Country, Year), .SDcols=c("supply", "export", "food", "feed", "import", "other", "processing", "prod", "dryMatter", "energy_kcal_pc", "stockvar")]
 	cb_fs_sugar[, c('Item', 'dryMatter', 'fcr', 'plumCropItem', 'plumDemandItem') := list("Sugar crops", dryMatter/2, NA, 'Sugar', 'Sugar')] # find mean of DM from sum	
 	cb_fs = rbind(cb_fs[!Item %in% c("Sugar cane", "Sugar beet"), list(Country, Year, Item, supply, export, food, feed, import,other, processing, prod, dryMatter, fcr, energy_kcal_pc, stockvar, plumCropItem, plumDemandItem)], cb_fs_sugar)
 	
-	# Get commodities produced from processing
-	cb_fs2 = merge (cb_fs[Item %in% itemMapping$resultant], itemMapping[, list(processed, Item=resultant, relPrice)], by="Item", all.x=TRUE, allow.cartesian=TRUE)	
-	cb_fs2[export>import, net_export:=export-import]
-	cb_fs2[is.na(net_export), net_export:=0]
-	cb_fs2[, used:=food+feed+other+net_export]
-	# cb_fs2[, relPrice:=1] THIS WOULD GIVE MASS ALLOCATION
-	cb_fs2[, weight:=relPrice*used/sum(relPrice*used), by=list(Country, Year, processed)]
-#	cb_fs2[import>export, energy_kcal_pc * prod / supply] # decrement energy to match domestic production, which is all that is being accounted for in produces from processing
-	
-	# Work out how commodities produced from processing are used
-	processWeights = cb_fs2[, list(food_rate = sum(food/used*weight), 
-					feed_rate = sum(feed/used*weight),
-					other_rate = sum(other/used*weight),
-					export_rate = sum(net_export/used*weight),
-					proc_kcal_pc = sum(energy_kcal_pc, na.rm=TRUE)), by=list(Country, Year, Item=processed)]
-	processWeights[is.na(food_rate), food_rate:=1]
-	
-	countries_missing_sugar_crops = cb_fs[!Country %in% cb_fs[Year == 2010 & Item == "Sugar crops"]$Country & Year == 2010 & Item == "Sugar (Raw Equivalent)"]$Country
-	sugar_energy = 0.00134996 #processedFull$cb_fs[Country == "World" & Year == 2010 & Item == "Sugar crops", list(supply /energy_kcal_pc)]/country_data[Country == "World" & Year == 2010]$population
-	sugar_missing = merge(processWeights[Country %in% countries_missing_sugar_crops & Item == "Sugar crops"], country_data[, list(Country, Year, population)], by=c("Country", "Year"))
-	sugar_missing[, supply:=proc_kcal_pc*sugar_energy*population]
-	sugar_missing = sugar_missing[, list(Country, Year, Item="Sugar crops", supply, export=supply*export_rate, food=supply*food_rate, feed=supply*feed_rate, import=0, other=supply*other_rate, 
-					processing=0, prod=supply, dryMatter=0.211, fcr=NA, energy_kcal_pc=0, stockvar=0, plumCropItem="Sugar", plumDemandItem="Sugar")]
-	cb_fs = rbind(cb_fs, sugar_missing)
-	
-	# Amend the FAO data based on the usage of commodities produced from processing
-	cb_fs3 = merge(cb_fs[!Item %in% itemMapping$resultant], processWeights, by=c("Country", "Year", "Item"), all.x=TRUE)
-	cb_fs3[is.na(energy_kcal_pc), energy_kcal_pc := 0]
-	cb_fs3[proc_kcal_pc > 0, energy_kcal_pc := energy_kcal_pc + proc_kcal_pc]
-	cb_fs3[proc_kcal_pc > 0, food := food + processing*food_rate]
-	cb_fs3[proc_kcal_pc > 0, feed := feed + processing*feed_rate]
-	cb_fs3[proc_kcal_pc > 0, other := other + processing*other_rate]
-	cb_fs3[proc_kcal_pc > 0, export := export + processing*export_rate]
-	cb_fs3[proc_kcal_pc > 0, processing:= 0]
+	itemMapping = fread("~/Documents/R_Workspace/landalloc/itemsProcessingMap.csv")
+	processedDt = merge(cb_fs[Item %in% itemMapping$resultant], itemMapping[, list(processed, Item=resultant)], by="Item", all.x=TRUE)	
+	extractionRates = extractionValues(cb_fs,processedDt)
+	itemMapping = merge(itemMapping, extractionRates[,list(processed,extractionRate)],by='processed' )
+	
+	#below converts processed goods into their raw equivalent values i.e. now columns for soybean oil are in soybean equiv. terms 
+	cb_fs3=merge(cb_fs,itemMapping[,list(Item=resultant,extractionRate)],by='Item',all=TRUE)
+	colsToAdj = c("supply", "export", "food", "feed", "import", "other", "processing", "stockvar", "dryMatter")
+	cb_fs3[!is.na(extractionRate), (colsToAdj) := lapply(.SD, function(x) x/extractionRate), by=list(Country, Year, Item), .SDcols=colsToAdj]
 	cb_fs3[, net_import := import-export]
-	cb_fs3[, c('food_rate', 'feed_rate', 'other_rate', 'export_rate', 'proc_kcal_pc', 'import', 'export') := NULL]
 	
-	histDt = getHistoricDataset(cb_fs3, country_data)[Year>=minYear]
+	#want to zero production of processed final goods as their production comes from the production of the raw good anyway 
+	cb_fs3[Item %in% itemMapping$resultant, prod := 0]
+	
+	histDt = getHistoricDataset(cb_fs3,country_data)[Year>=minYear]
 	ssp=getSsp(country_data)
 	com_curves = getCommodityCurves(histDt)
 	baseConsumption = getBaseConsumption(year,histDt,ssp)
@@ -248,4 +238,4 @@ processData = function(minYear=1960, year=2010){
 					projectedConsumption=projectedConsumption, cb_fs=cb_fs3, livestockNumbers=livestockNumbers))
 }
 
-processedFull = processData()
\ No newline at end of file
+processedFull = processData()
diff --git a/PLUM/createShockFiles.R b/PLUM/createShockFiles.R
index a831a3411b3ba905b96bb94adefac5e2c2901b5c..4506dd9d5638a0b486530bb2478426b8d14602f5 100644
--- a/PLUM/createShockFiles.R
+++ b/PLUM/createShockFiles.R
@@ -11,7 +11,7 @@ createShockFile = function(scenario){
 	
 	paramValues = data.table(Year = seq(2010,2100,1),
 			trendYearValue = c(rep(0,16),seq(0.2,1,0.2),rep(1,70)),
-			TRADE_BARRIERS_MULTIPLIER =scenario[,TRADE_BARRIERS_MULTIPLIER], 
+			TRADE_BARRIER_MULTIPLIER =scenario[,TRADE_BARRIER_MULTIPLIER], 
 			TRANSPORT_COST = scenario[,TRANSPORT_COST],
 			TECHNOLOGY_CHANGE_ANNUAL_RATE = scenario[,TECHNOLOGY_CHANGE_ANNUAL_RATE],
 			MEAT_EFFICIENCY = scenario[,MEAT_EFFICIENCY],
@@ -27,12 +27,12 @@ createShockFile = function(scenario){
 	paramValues[, names(paramValues) := lapply(.SD, as.numeric)]
 	
 	shocksHappen = data.table(SSP_Scenario = scenario[,SSP_SCENARIO], Scenario = scenario[,Scenario], Year = seq(2010,2100,1), 
-			climate =as.logical(rbinom(91,size=1,prob = shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),climate])), 
+			climate =as.logical(rbinom(91,size=1,prob = climateShockProbs[rcp == scenario[,RCP],climate])), 
 			cyber=as.logical(rbinom(91,size=1,prob = shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),cyber])), 
 			financial =as.logical(rbinom(91,size=1,prob = shockProbs[ssp == substr(scenario[,SSP_SCENARIO],1,4),financial])))
 	
 	shockValues = data.table(Year = rep(seq(2010,2100,1),6),shock = rep(c("protectionism","automation","diet","cyber","climate","financial"), each=91), 
-			TRADE_BARRIERS_MULTIPLIER = 0, 
+			TRADE_BARRIER_MULTIPLIER = 0, 
 			TRANSPORT_COST = 0,
 			TECHNOLOGY_CHANGE_ANNUAL_RATE = 0,
 			MEAT_EFFICIENCY = 0,
@@ -85,6 +85,20 @@ createShockFile = function(scenario){
 		}
 	}
 	
+	
+	yieldshocks = NULL
+	for(year in shocksHappen[,Year]){
+		if(shocksHappen[Year == year, climate]){
+			
+			mapFilename = sample(climateMaps, 1)
+			shockFactor = runif(1,0.1,0.5)
+			yieldshocks = rbind(yieldshocks, cbind(year,mapFilename,crop='maize',shockFactor))
+			yieldshocks = rbind(yieldshocks, cbind(year,mapFilename,crop='rice',shockFactor))
+			yieldshocks = rbind(yieldshocks, cbind(year,mapFilename,crop='cereal',shockFactor))
+		}		
+		
+	}
+	
 	#if negative set to zero, or need to change shock distributions so zero can't happen 
 	for(col in names(paramValues)) set(paramValues, i=which(paramValues[[col]]<0), j=col, value=0)
 	
@@ -93,7 +107,8 @@ createShockFile = function(scenario){
 	if (dir.exists(file.path(baseOutputDir,ensDir))) {
 		print(paste(ensDir, "exists.  Creating shock file"))
 		write.csv(paramValues, file.path(baseOutputDir,ensDir,"shocks.csv"), row.names=FALSE, quote=FALSE)
-		write.csv(shockValues, file.path(baseOutputDir,ensDir,"shocksRecord.csv"), row.names=FALSE, quote=FALSE)
+		write.csv(shockValues, file.path(baseOutputDir,ensDir,paste0("shocksRecord_",scenario[,Scenario],".csv")), row.names=FALSE, quote=FALSE)
+		write.csv(yieldshocks, file.path(baseOutputDir,ensDir,paste0("yieldshocks.csv")), row.names=FALSE, quote=FALSE)
 	}
 	else{
 		print(paste(ensDir, "does not exist.  Stopping"))
@@ -112,15 +127,19 @@ drawValue <- function(shockType, paramName){
 	randomValue
 }
 
+plumDirectory = "/exports/csce/eddie/geos/groups/LURG/models/PLUM"
 baseDataDir = "/exports/csce/eddie/geos/groups/LURG/models/PLUM/plumv2/data"
 baseOutputDir="/exports/csce/eddie/geos/groups/LURG/models/PLUM/output" # "~/Downloads"
+climateMaps = basename(list.files(path = file.path(baseDataDir,"yieldshockmaps"), full.names = TRUE, recursive = TRUE))
 
 shockDistFile = fread(file.path(baseDataDir, "shockDistributions.csv"))
 shockProbs = fread(file.path(baseDataDir,"shockProbabilities.csv"))
+climateShockProbs = fread(file.path(baseDataDir,"climateShockProbabilities.csv"))
+
 
 ensembleFile = commandArgs(trailingOnly = TRUE)[1]
 
-scenarioTable = fread(file.path(baseDataDir,"sims",ensembleFile))
+scenarioTable = fread(file.path(plumDirectory,ensembleFile))
 
 for(i in 1:nrow(scenarioTable)){
 #first if column exists
diff --git a/PLUM/debugging.Rmd b/PLUM/debugging.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..f5afb18b8eae5a6a889979e3c3d1d9684bba1fed
--- /dev/null
+++ b/PLUM/debugging.Rmd
@@ -0,0 +1,338 @@
+---
+title: "PLUM debugging script"
+author: "Roslyn Henry"
+date: "24 April 2019"
+output:
+  html_document:
+    fig_width: 6
+    fig_height: 4
+---
+
+```{r setup, include=FALSE}
+
+library(plotly)
+library(ggplot2)
+library(gridExtra)
+library(rmarkdown)
+library(data.table)
+library(raster)
+
+require(reshape2)
+require(maptools)
+
+
+plum_data_dir = "C:\\Users\\rhenry2\\workspace\\PLUM_src\\plumv2\\plumv2\\data"
+Directory='C:\\Users\\rhenry2\\Desktop\\s_newds4_rh'
+
+shpFile = crop(readShapePoly("C:\\Users\\rhenry2\\Desktop\\Biodiversity_data\\NaturalEarthShapefile\\NaturalEarthShapefile\\ne_110m_admin_0_countries_lakes.shp", proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")), extent(-180,180, -55, 90))
+
+```
+
+
+```{r functions, include=FALSE}
+
+
+plotLandUse=function(directory){
+	lc=fread(paste0(directory,"\\lc.txt"))
+	
+	h_melt = melt(lc[, list(Year,natural=Natural+ManForest+UnmanForest, cropland=Cropland, pasture=Pasture, fert=FertCrop, irrig=IrrigCrop)], id.vars=c('Year'), variable.name='variable', value.name="Value")
+	
+	cRange = range(h_melt[variable=="cropland"]$Value)
+	pRange = range(h_melt[variable=="pasture"]$Value)
+	fRange = range(h_melt[variable=="fert"]$Value)
+	wRange = range(h_melt[variable=="irrig"]$Value)
+	nRange = range(h_melt[variable=="natural"]$Value)
+	
+	h_melt[variable=="cropland", variable:="a) Cropland area\n(Mha)"]
+	h_melt[variable=="pasture", variable:="b) Pasture area\n(Mha)"]
+	h_melt[variable=="natural", variable:="c) Natural area\n(Mha)"]
+	h_melt[variable=="fert", variable:="d) Nitrogen applied\n(Mt)"]
+	h_melt[variable=="irrig", variable:="e) Irrigation water\nwithdrawn (km\u00b3)"]
+	h_melt[, variable:=factor(variable, levels=c("a) Cropland area\n(Mha)","b) Pasture area\n(Mha)", "c) Natural area\n(Mha)", "d) Nitrogen applied\n(Mt)", "e) Irrigation water\nwithdrawn (km\u00b3)"))]
+	
+	h_avg = h_melt[, list(sdev = sd(Value), avg = median(Value)), by=list(Year, variable)]
+	
+	getConePlotForType = function (dt, dType, scaleRange, colours=c("skyblue4","grey","green")) {
+		d = h_avg[variable==dType]
+		
+		p1 = ggplot(d) +
+				geom_line(aes(x=Year, y=avg)) +
+				labs(y=dType) +
+				scale_y_continuous(limits=scaleRange) +
+				scale_x_continuous(minor_breaks = seq(2010, 2100, 10), breaks = seq(2010, 2100, 20)) + 
+				theme_minimal() +geom_ribbon(aes(x=Year, ymin=avg-sdev,ymax=avg+sdev), alpha=0.05)
+		
+		p1
+	}
+	
+	
+	cp= getConePlotForType(h_avg, dType="a) Cropland area\n(Mha)", scaleRange=cRange)
+	pp = getConePlotForType(h_avg, dType="b) Pasture area\n(Mha)", scaleRange=pRange)
+	np = getConePlotForType(h_avg, dType="c) Natural area\n(Mha)", scaleRange=nRange)
+	fp = getConePlotForType(h_avg, dType="d) Nitrogen applied\n(Mt)", scaleRange=fRange)
+	wp = getConePlotForType(h_avg, dType="e) Irrigation water\nwithdrawn (km\u00b3)", scaleRange=wRange)
+	
+	grid.arrange(cp, pp, np , fp, wp, ncol=2, widths=c(1, 1), padding=0)
+}
+
+getFeedComparison = function(directory){
+	
+	fbs = fread(paste0(directory,"\\fbs.txt"))
+	fbs[Crop == 'pulses', Crop := 'Pulses']
+	fbs[Crop == 'oilcrops', Crop := 'Oilcrops']
+	fbs[Crop == 'starchyRoots', Crop := 'Starchy Roots']
+	fbs[Crop == 'sugar', Crop := 'Sugar']
+	fbs[Crop == 'fruitveg', Crop := 'FruitVeg']
+	fbs[Crop == 'maize', Crop := 'MaizeMilletSorghum']
+	fbs[Crop == 'wheat', Crop := 'WheatBarleyOats']
+	fbs[Crop == '"rice', Crop := 'Rice (Paddy Equivalent)']
+	
+	fao = processedFull$cb_fs[is.na(fcr) & Country=="World" & Year == 2010, list(faoFeed=sum(feed), faoFeedDM=sum(feed*dryMatter)), by=list(Crop=plumCropItem)]
+	plum=merge(fbs, itemDetails[, list(dryMatter=mean(dryMatter)), by=list(Crop=plumCropItem)], by="Crop")[Year == 2010, list(mono=sum(MonogastricsFeed), rum=sum(RuminantsFeed), monoDM=sum(MonogastricsFeed * dryMatter), rumDM=sum(RuminantsFeed * dryMatter)), by=Crop]
+	feedComparison=merge(fao, plum, by="Crop")
+	feedComparison[, diff:=mono+rum-faoFeed]
+	feedComparison[, diffDM:=monoDM+rumDM-faoFeedDM]
+	feedComparison[, list(faoFeed=sum(faoFeed), faoFeedDM =sum(faoFeedDM))]
+	feedComparison
+
+}
+
+
+plotPriceFileOutput=function(dt, measure, crops, ylabel){
+	
+		ggplot(dt[type==measure & Crop %in% crops] , aes(x=Year, y=Index_value)) + geom_line(size=1)+
+				facet_grid(~Crop, scales="free_y") + scale_y_continuous() + ylab(ylabel)+
+				theme_minimal()+ scale_x_continuous() +
+				theme(legend.position="bottom",legend.title=element_blank(),
+						plot.margin=margin(1,1,1,1, 'in'),
+						panel.border = element_rect(colour = "grey10", fill=NA, size=0.3),
+						axis.text=element_text(size=20),
+						axis.title.x=element_text(size=20,face="bold"),
+						axis.title.y=element_text(size=20, face="bold"),
+						legend.text=element_text(size=20),
+						strip.text.x = element_text(size =20),
+						strip.text.y = element_blank())
+
+}
+
+rasterFromFile = function(lcType, fileName, offset=0.0) {
+	plumDt=data.table(read.table(fileName, header=TRUE))
+	plumDt[, c('Lon', 'Lat') := list(Lon-offset, Lat-offset)]
+	rasterFromXYZ(plumDt[, list(Lon, Lat, get(lcType))], crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
+}
+
+plotDiff = function(r1, r2, label) {
+	
+	e = intersect(extent(r1), extent(r2))
+	rDiff = crop(r2, e)-crop(r1, e)
+	
+	breakpoints <- c(seq(from=-1,to=1,by=0.02))#0,seq(from=0,to=0.1,by=0.02))
+	colfuncR <- colorRampPalette(c("firebrick1","gray97","blue"))
+	colors <- c(colfuncR(100))
+	
+	plot(rDiff,
+			col=colors,breaks=breakpoints, box=FALSE,axes = FALSE,
+			horizontal=TRUE,legend.mar=3,legend.lab=label,
+			axis.args=list(at=seq(-1, 1, 0.2), labels=seq(-1, 1, 0.2), 
+					cex.axis=1,  mgp=c(1, .1, 0)))
+	abline(h=23.5, col="black", lty='dotted')
+	abline(h=-23.5, col="black", lty='dotted')	
+	plot(shpFile, add=TRUE,border="black", lwd=0.05)
+	#mtext(label, side=3, line=-23, cex=1.75)
+	
+}
+
+plotCandPChanges = function(directory, startYr, endYr){
+
+	par(mfrow=c(2, 1), mar=c(2, 2, 2, 2), oma=c(3,3,3,3))
+	cs = rasterFromFile("CROPLAND", fileName=file.path(directory, startYr, "LandCoverFract.txt"))
+	cf = rasterFromFile("CROPLAND", fileName=file.path(directory, endYr, "LandCoverFract.txt"))
+	ps = rasterFromFile("PASTURE", fileName=file.path(directory, startYr, "LandCoverFract.txt"))
+	pf = rasterFromFile("PASTURE", fileName=file.path(directory, endYr, "LandCoverFract.txt"))
+	
+	plotDiff(cs, cf, 'Cropland change fraction')
+	plotDiff(ps, pf, 'Pasture change fraction')
+	
+}
+
+```
+
+
+```{r prepare_data, include=FALSE}
+##format plum prices file for plotting below 
+prices = fread(paste0(Directory,'\\prices.txt'))
+names(prices)= gsub(" ", "_", names(prices))
+
+price_melt = melt(prices, measure.vars=c("Imports_(Mt)","Exports_(Mt)",
+				"New_export_price","Stock_Levels_(Mt)"), variable.name='type', value.name="Index_value")
+
+##format plum demand system output file for plotting below 
+countryDemandOpt = fread(paste0(Directory,'\\countryDemandOpt.txt'))
+countryDemandOpt[commodity == 'cereals', Item := 'Cereals']
+countryDemandOpt[commodity == 'starchyRoots', Item := 'Starchy Roots']
+countryDemandOpt[commodity == 'oilcropspulses', Item := 'OilcropsPulses']
+countryDemandOpt[commodity == 'sugar', Item := 'Sugar']
+countryDemandOpt[commodity == 'fruitveg', Item := 'FruitVeg']
+countryDemandOpt[commodity == 'ruminants', Item := 'Ruminants']
+countryDemandOpt[commodity == 'monogastrics', Item := 'Monogastrics']
+
+countryDemandOpt1 = countryDemandOpt[commodity != 'nonfood',list(Item = 'All', plumRebasedMt=sum(plumRebased), plumNotRebaseMt=sum(plumNotRebase),rebasedKcal=sum(rebasedKcal),notRebasedKcal=sum((subsistence+discretionary)*2000)), by=list(country,year)]
+countryDemandOpt = rbind(countryDemandOpt[,list(Country=country, Year=year, plumDemandItem=Item,plumRebasedMt=plumRebased,plumNotRebaseMt=plumNotRebase, rebasedKcal=rebasedKcal, notRebasedKcal=(subsistence+discretionary)*2000)],
+		countryDemandOpt1[,list(Country=country, Year=year, plumDemandItem=Item,plumRebasedMt,plumNotRebaseMt, rebasedKcal, notRebasedKcal)])
+
+#processedFull comes from readFAO script. Need also to read in commonFunctions script
+processedFull = processData()
+countryGroups =fread(paste0(plum_data_dir,"\\country_groups.csv"), header=TRUE)
+
+caloriesPerCountryFAO = processedFull$cb_fs[Year==2010 & Country %in% countryGroups[,Country], list(caloriesPc=sum(energy_kcal_pc, na.rm=TRUE), food=sum(food,na.rm=TRUE)), by=c('Country','plumDemandItem') ]
+caloriesPerCountryFAO = rbind(caloriesPerCountryFAO,caloriesPerCountryFAO[,list(plumDemandItem='All', caloriesPc=sum(caloriesPc), food=sum(food)), by='Country'])
+cals = merge(processedFull$cb_fs[Year == 2010],processedFull$country_data[Year == 2010, list(Country,population, gdp_pc)], by=c('Country'))
+cals[!is.na(fcr), c('food') := list(food * fcr)]
+cals2 = cals[,list(caloriesPc=sum(energy_kcal_pc, na.rm=TRUE), baseCpc=sum(food/population,na.rm=TRUE)), by=c('Country','plumDemandItem') ]
+
+cals3=merge(cals2,countryDemandOpt[Year==2010,list(Country, plumDemandItem,notRebasedKcal,rebasedKcal, plumRebasedMt, plumNotRebaseMt)], by=c('Country','plumDemandItem'))
+
+cals3[,percKcal := notRebasedKcal/caloriesPc]
+cals3[,percMt := plumNotRebaseMt/baseCpc]
+cals3[,percMtRebased := plumRebasedMt/baseCpc]
+cals3[,percKcalRebased := rebasedKcal/caloriesPc]
+
+
+```
+
+**Global land use change over time.**
+
+```{r land_use_change, include=TRUE, echo=FALSE, fig.width=12, fig.height=8}
+	plotLandUse(directory=Directory)
+```
+
+**Land use change over between 2010 and 2100.**
+
+```{r spatial_land_use_change, include=TRUE, echo=FALSE, fig.width=12, fig.height=8}
+	plotCandPChanges(directory=Directory, 2010,2100)
+```
+
+**Comparing FAO feed use in 2010 to PLUM feed use in 2010.** 
+
+
+```{r  feed_use, include=TRUE, echo=FALSE}
+	feedComparison=getFeedComparison(directory=Directory)
+	feedComparison
+	feedComparison[, list(diffSum=sum(diff), diffSumDM=sum(diffDM))]
+```
+
+**Global price over time.** 
+
+```{r price_over_time, include=TRUE, echo=FALSE, fig.width=12, fig.height=8}
+
+plotPriceFileOutput(dt=price_melt, measure='New_export_price', crops=c("ruminants","monogastrics","fruitveg"), 
+		ylabel='Global price')
+plotPriceFileOutput(dt=price_melt, measure='New_export_price', crops=c("wheat","maize","rice"), 
+		ylabel='Global price')
+plotPriceFileOutput(dt=price_melt, measure='New_export_price', crops=c("oilcrops","pulses","starchyRoots","sugar"), 
+		ylabel='Global price')
+
+```
+
+**Global imports over time.**
+
+
+```{r imports_over_time, include=TRUE, echo=FALSE, fig.width=12, fig.height=8}
+
+plotPriceFileOutput(dt=price_melt, measure='Imports_(Mt)', crops=c("ruminants","monogastrics","fruitveg"), 
+		ylabel='Imports (Mt)')
+plotPriceFileOutput(dt=price_melt, measure='Imports_(Mt)', crops=c("wheat","maize","rice"), 
+		ylabel='Imports (Mt)')
+plotPriceFileOutput(dt=price_melt, measure='Imports_(Mt)', crops=c("oilcrops","pulses","starchyRoots","sugar"), 
+		ylabel='Imports (Mt)')
+
+```
+
+**Calorie intake over time.**
+
+```{r Cals_rebased_over_time, include=TRUE, echo=FALSE, fig.width=16, fig.height=8}
+
+ggplotly(ggplot(countryDemandOpt[plumDemandItem %in% c('All','Ruminants','Monogastrics')],
+aes(x=Year,y=rebasedKcal, col=Country)) + geom_line()+
+		facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + ylab("Calories consumed")+ xlab("Year")+
+		theme(axis.title.y = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0)),
+				axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))))
+
+ggplotly(ggplot(countryDemandOpt[plumDemandItem %in% c('OilcropsPulses','Sugar','FruitVeg')],
+						aes(x=Year,y=rebasedKcal, col=Country)) + geom_line()+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + ylab("Calories consumed")+ xlab("Year")+
+				theme(axis.title.y = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0)),
+						axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))))
+
+ggplotly(ggplot(countryDemandOpt[plumDemandItem %in% c('Starchy Roots','Cereals')],
+						aes(x=Year,y=rebasedKcal, col=Country)) + geom_line()+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + ylab("Calories consumed")+ xlab("Year")+
+				theme(axis.title.y = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0)),
+						axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))))
+
+```
+
+**Calorie intake projected by demand system over time.**
+
+```{r Cals_ds_over_time, include=TRUE, echo=FALSE, fig.width=16, fig.height=8}
+
+
+ggplotly(ggplot(countryDemandOpt[plumDemandItem %in% c('All','Ruminants','Monogastrics')],
+						aes(x=Year,y=notRebasedKcal, col=Country)) + geom_line()+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + ylab("Calories consumed demand system")+ xlab("Year")+
+				theme(axis.title.y = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0)),
+						axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))))
+
+ggplotly(ggplot(countryDemandOpt[plumDemandItem %in% c('OilcropsPulses','Sugar','FruitVeg')],
+						aes(x=Year,y=notRebasedKcal, col=Country)) + geom_line()+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + ylab("Calories consumed demand system")+ xlab("Year")+
+				theme(axis.title.y = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0)),
+						axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))))
+
+ggplotly(ggplot(countryDemandOpt[plumDemandItem %in% c('Starchy Roots','Cereals')],
+						aes(x=Year,y=notRebasedKcal, col=Country)) + geom_line()+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + ylab("Calories consumed demand system")+ xlab("Year")+
+				theme(axis.title.y = element_text(margin = margin(t = 0, r = 0, b = 20, l = 0)),
+						axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l = 0))))
+```
+
+
+**Calorie intake in 2010 projected by demand system compared to FAO food
+supply data in 2010.**
+
+```{r Cals_ds_vs_FAO, include=TRUE, echo=FALSE, fig.width=16, fig.height=8}
+
+
+ggplotly(ggplot(cals3[plumDemandItem %in% c('All','Ruminants','Monogastrics')], aes(x=caloriesPc,y=notRebasedKcal, col=Country)) + geom_point(size=2)+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + scale_x_continuous()
+				+ xlab("FAO kcal per capita")+
+				ylab("Demand system kcal per capita"))
+
+ggplotly(ggplot(cals3[plumDemandItem %in% c('OilcropsPulses','Sugar','FruitVeg')], aes(x=caloriesPc,y=notRebasedKcal, col=Country)) + geom_point(size=2)+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + scale_x_continuous()
+				+ xlab("FAO kcal per capita")+
+				ylab("Demand system kcal per capita"))
+
+ggplotly(ggplot(cals3[plumDemandItem %in% c('Starchy Roots','Cereals')], aes(x=caloriesPc,y=notRebasedKcal, col=Country)) + geom_point(size=2)+
+				facet_grid(.~plumDemandItem, scales="free_y") + scale_y_continuous() + scale_x_continuous()
+				+ xlab("FAO kcal per capita")+
+				ylab("Demand system kcal per capita"))
+
+```
+
+**Looking at variation against FAO data with rebased and original demand
+output**
+
+```{r variation_calories, include=TRUE, echo=FALSE}
+
+cals3[,list(mPercKcal = mean(percKcal, na.rm=TRUE), 
+				sdPercKcal = sd(percKcal, na.rm=TRUE), 
+				mPercKcalRebased = mean(percKcalRebased, na.rm=TRUE), 
+				sdPercKcalRebased = sd(percKcalRebased, na.rm=TRUE), 
+				mPercMt = mean(percMt, na.rm=TRUE), 
+				sdPercMt=sd(percMt, na.rm=TRUE),
+				mpercMtRebased = mean(percMtRebased, na.rm=TRUE), 
+				sdpercMtRebased=sd(percMtRebased, na.rm=TRUE)), 
+		by='plumDemandItem']
+```
+
diff --git a/itemsInputs.csv b/itemsInputs.csv
index 46ffd620f530a5ebcb7c9a8d4adc64d583f0deb6..d093550528744e0906efc8b68bf9f0abcd64f8b0 100644
--- a/itemsInputs.csv
+++ b/itemsInputs.csv
@@ -35,11 +35,11 @@ Maize Germ Oil,FALSE,FALSE,1,Cereals,MaizeMilletSorghum,,Cereals
 Meat Other,FALSE,TRUE,0.33,Meat,Ruminants,15,Ruminants
 Milk - Excluding Butter,FALSE,TRUE,0.119,Milk,Ruminants,0.7,Ruminants
 Millet,TRUE,FALSE,0.902,Cereals,MaizeMilletSorghum,,Cereals
-Molasses,FALSE,FALSE,0.73,Sugar and sweeteners,FruitVeg,,FruitVeg
+Molasses,FALSE,FALSE,0.73,Sugar and sweeteners,Sugar,,Sugar
 Mutton & Goat Meat,FALSE,TRUE,0.405,Meat,Ruminants,15,Ruminants
 Nuts,TRUE,FALSE,0.9,Oilseeds and pulses,Oilcrops,,OilcropsPulses
 Oats,TRUE,FALSE,0.879,Cereals,WheatBarleyOats,,Cereals
-Oil palm fruit,TRUE,FALSE,0.34,Oilseeds and pulses,Oilcrops,,OilcropsPulses
+Palm Oil,TRUE,FALSE,1,Oilseeds and pulses,Oilcrops,,OilcropsPulses
 Oilcrops Oil Other,FALSE,FALSE,1,Oilseeds and pulses,Oilcrops,,OilcropsPulses
 Oilcrops Other,TRUE,FALSE,0.88,Oilseeds and pulses,Oilcrops,,OilcropsPulses
 Oilseed Cakes Other,FALSE,FALSE,0.9,Oilseeds and pulses,Oilcrops,,OilcropsPulses
diff --git a/itemsProcessingMap.csv b/itemsProcessingMap.csv
index d711ded88c5fc7f6f17192a1b01577e350e7ec34..7ba004954efd4982b0d0de0f3d92f5a227f90d2f 100644
--- a/itemsProcessingMap.csv
+++ b/itemsProcessingMap.csv
@@ -13,8 +13,8 @@ Maize,Sweeteners Other,2
 Oilcrops Other,Oilcrops Oil Other,1
 Oilcrops Other,Oilseed Cakes Other,3
 Olives (including preserved),Olive Oil,1
-Palm Kernels,Palmkernel Cake,1
-Palm Kernels,Palmkernel Oil,3
+Palm kernels,Palmkernel Oil,3
+Palm kernels,Palmkernel Cake,1
 Potatoes,Beverages Fermented,1
 Rape and Mustardseed,Rape and Mustard Cake,1
 Rape and Mustardseed,Rape and Mustard Oil,3