Prepare global soil carbon point data set

ISRICWorldSoil · Oct 11, 2017 · 5ae9930 · 5ae9930
1 parent 4b89ec0
commit 5ae9930
Show file tree

Hide file tree

Showing 5 changed files with 93 additions and 6 deletions.
diff --git a/grids/aggregated/preview_SoilGrids.R b/grids/aggregated/preview_SoilGrids.R
@@ -96,6 +96,12 @@ write.csv(Tocs30.sum, "Status_OCS_30cm_per_Country.csv")
 ## Total stock 0-30cm:
 sum(Tocs30.sum$OCS_mT[!is.na(Tocs30.sum$Value)])
 
+## Compare with 1km data:
+grid1km.sin = readGDAL("/data/aggregated/1km/OCSTHA_M_30cm_1km_sin.tif")
+sum(grid1km.sin$band1*1000^2/10000, na.rm=TRUE)/1e6
+## 1360 Pg C
+rm(grid1km.sin); gc()
+
 ## List of property maps:
 tif.lst <- list.files(path="/data/GEOG", pattern=glob2rx("*.tif$"), full.names=TRUE, recursive=TRUE)
 ## 318

diff --git a/grids/models/SPROPS/ORC_summary_per_TAXOUSDA.csv b/grids/models/SPROPS/ORC_summary_per_TAXOUSDA.csv
@@ -0,0 +1,66 @@
+"","TAXOUSDA","M_ORC_30","M_ORC_100","M_BLD_30","M_BLD_100","M_OCD_30","N_horizons","TAXOUSDA_name","AREA"
+"1","5",423,436,246,244,86,4694,"Histels",1850412
+"2","6",186,135,678,729,73,1938,"Turbels",3315040
+"3","7",130,87,743,930,64,2240,"Orthels",8716873
+"4","10",296,266,647,721,161,33,"Folists",88
+"5","11",343,325,331,348,72,916,"Fibrists",389985
+"6","12",307,334,326,315,87,2140,"Hemists",113966
+"7","13",252,240,592,590,150,1853,"Saprists",61591
+"8","15",46,29,866,1065,67,5367,"Aquods",1020821
+"9","16",127,95,840,922,89,691,"Cryods",33575
+"10","17",92,59,733,915,69,259,"Humods",2505
+"11","18",80,55,779,1032,53,22639,"Orthods",4689232
+"12","20",123,103,890,1046,87,160,"Aquands",311
+"13","21",104,62,772,938,55,1420,"Cryands",97387
+"14","22",16,11,1430,1435,19,18,"Torrands",9
+"15","23",130,76,760,965,61,1389,"Xerands",32040
+"16","24",65,35,886,1148,45,1061,"Vitrands",35743
+"17","25",112,81,856,932,75,254,"Ustands",1173
+"18","26",133,87,928,1057,78,2718,"Udands",456831
+"19","30",34,26,1197,1274,35,90,"Aquox",167
+"20","31",20,13,1295,1325,29,1238,"Torrox",307406
+"21","32",20,16,1210,1217,21,30215,"Ustox",11049085
+"22","33",57,45,1040,1075,57,122,"Perox",1901
+"23","34",25,19,1079,1160,35,10239,"Udox",4514812
+"24","40",22,15,1507,1575,33,1328,"Aquerts",64963
+"25","42",11,8,1646,1698,17,671,"Xererts",29374
+"26","43",12,9,1563,1631,22,325,"Torrerts",3166
+"27","44",13,10,1506,1531,18,8482,"Usterts",3035746
+"28","45",20,13,1706,1762,32,822,"Uderts",25311
+"29","50",23,16,1257,1320,18,45,"Cryids",227
+"30","51",6,4,1404,1454,6,150,"Salids",4028
+"31","52",7,5,1459,1482,8,437,"Durids",35176
+"32","53",7,4,1307,1308,7,13215,"Gypsids",12083
+"33","54",7,5,1507,1532,9,10335,"Argids",2314899
+"34","55",9,7,1452,1478,11,3897,"Calcids",964779
+"35","56",7,5,1426,1435,12,2810,"Cambids",524287
+"36","60",23,14,1360,1504,38,1557,"Aquults",301237
+"37","61",45,33,993,1154,46,5308,"Humults",564160
+"38","62",23,16,1290,1387,27,51132,"Udults",16477412
+"39","63",14,12,1334,1341,17,7196,"Ustults",1695001
+"40","64",56,32,1080,1235,46,288,"Xerults",1086
+"41","69",33,22,1237,1356,32,10529,"Borolls",1337910
+"42","70",53,32,1059,1275,35,1365,"Albolls",60588
+"43","71",38,23,1213,1388,41,12598,"Aquolls",2980442
+"44","72",39,30,965,1133,43,251,"Rendolls",14735
+"45","73",25,19,1301,1392,30,12574,"Xerolls",2381368
+"46","74",43,25,1213,1349,35,1650,"Cryolls",449008
+"47","75",20,14,1386,1441,22,34646,"Ustolls",5367823
+"48","76",20,13,1426,1481,27,25004,"Udolls",2109178
+"49","80",16,10,1454,1562,22,15172,"Aqualfs",422040
+"50","81",53,28,1170,1355,39,525,"Cryalfs",50682
+"51","82",12,9,1425,1463,14,34414,"Ustalfs",5223450
+"52","83",34,23,1288,1438,32,6016,"Xeralfs",406989
+"53","84",18,12,1338,1472,27,84597,"Udalfs",3660278
+"54","85",77,48,936,1183,63,9315,"Udepts",907814
+"55","89",72,50,889,1114,38,25984,"Ochrepts",6977485
+"56","90",68,45,933,1131,51,10542,"Aquepts",898137
+"57","92",71,42,922,1123,47,1032,"Cryepts",137880
+"58","93",21,13,1414,1485,20,1355,"Ustepts",124680
+"59","94",58,41,874,1062,65,2178,"Xerepts",286309
+"60","95",37,28,821,946,37,11725,"Aquents",743305
+"61","96",19,14,1156,1198,27,3263,"Arents",238492
+"62","97",12,9,1404,1420,9,22249,"Psamments",18041699
+"63","98",22,15,1089,1176,25,10900,"Fluvents",1352518
+"64","99",22,17,1262,1311,24,22118,"Orthents",12301141
+"65",NA,40,26,1156,1277,34,123382,NA,NA
diff --git a/grids/models/SPROPS/m.SPROPS.R b/grids/models/SPROPS/m.SPROPS.R
@@ -29,7 +29,7 @@ library(plotKML)
 library(R.utils)
 library(GSIF)
 library(parallel)
-library(doParallel)
+#library(doParallel)
 
 plotKML.env(convert="convert", show.env=FALSE)
 gdalwarp = "gdalwarp"
@@ -91,6 +91,7 @@ hexbinplot(ovA$BLD~ovA$ORCDRC, colramp=colorRampPalette(pal), xlab="Organic carb
 ## Fill in gaps in soil BLD - we can make a PedoTransfer function (which includes soil type predictions)
 ov.tax = raster::extract(raster("/data/GEOG/TAXOUSDA_250m_ll.tif"), SPROPS.pnts)
 dfsB = plyr::join(ovA[,c("SOURCEID", "BLD", "ORCDRC", "CLYPPT", "SNDPPT", "PHIHOX", "DEPTH.f")], data.frame(SOURCEID=SPROPS.pnts$SOURCEID, TAXOUSDA=as.factor(ov.tax)))
+## data to fit BLD PTF:
 sel.ovA = complete.cases(dfsB[,c("ORCDRC", "BLD", "CLYPPT", "SNDPPT", "PHIHOX", "DEPTH.f", "TAXOUSDA")])
 summary(sel.ovA)
 dfsB = dfsB[sel.ovA,]
@@ -175,7 +176,20 @@ View(ovA[which(ovA$OCDENS>200),c("SOURCEID","BLD","ORCDRC","BLD.f","OCDENS")])
 hexbinplot(log1p(ovA$OCDENS)~log1p(ovA$ORCDRC), colramp=colorRampPalette(pal), xlab="log - Organic carbon (permille)", ylab="log - Organic carbon density (kg/cubic-m)", type="g", lwd=1, lcex=8, inner=.2, cex.labels=.8, xbins=30, ybins=30, panel=pfun, colorcut=c(0,0.01,0.03,0.07,0.15,0.25,0.5,0.75,1))
 ## treshold at ca 12%
 expm1(4.8)
-
+## Average ORC, BLD per soil type:
+s.ovA = plyr::join(ovA[,c("SOURCEID", "SOURCEDB", "TIMESTRR", "LONWGS84", "LATWGS84", "HZDTXT", "UHDICM", "LHDICM", "BLD", "BLD.f", "ORCDRC", "OCDENS", "CLYPPT", "SNDPPT", "PHIHOX", "DEPTH.f")], data.frame(SOURCEID=SPROPS.pnts$SOURCEID, TAXOUSDA=as.factor(ov.tax)))
+sum.dfsB = plyr::ddply(s.ovA, .(TAXOUSDA), summarize, M_ORC_30=round(mean(ORCDRC[DEPTH.f<30], na.rm=TRUE)), M_ORC_100=round(mean(ORCDRC[DEPTH.f<100], na.rm=TRUE)), M_BLD_30=round(mean(BLD.f[DEPTH.f<30], na.rm=TRUE)), M_BLD_100=round(mean(BLD.f[DEPTH.f<100], na.rm=TRUE)), M_OCD_30=round(mean(OCDENS[DEPTH.f<30], na.rm=TRUE)), N_horizons=sum(!is.na(ORCDRC)))
+TAXOUSDA.leg <- read.csv("../TAXOUSDA/TAXOUSDA_legend.csv")
+sum.dfsB$TAXOUSDA_name = join(data.frame(Number=sum.dfsB$TAXOUSDA), TAXOUSDA.leg)$Group
+## Total areas per soil type:
+grid1km.sin = readGDAL("/data/aggregated/1km/TAXOUSDA_1km_sin.tif")
+## 20015 rows and 40030 columns
+sum.USDA = summary(as.factor(grid1km.sin$band1))
+sum.dfsB$AREA = join(sum.dfsB, data.frame(TAXOUSDA=attr(sum.USDA, "names"), AREA=sum.USDA))$AREA
+write.csv(sum.dfsB, "ORC_summary_per_TAXOUSDA.csv")
+saveRDS(s.ovA, "global_soil_carbon_points.rds")
+
+## Save a copy of all soil point data / regression matrix:
 write.csv(ovA, file="ov.SPROPS_SoilGrids250m.csv")
 unlink("ov.SPROPS_SoilGrids250m.csv.gz")
 gzip("ov.SPROPS_SoilGrids250m.csv")
@@ -367,7 +381,7 @@ check_readGDAL = function(x){
   }
   return(out)
 }
-## takes 120 mins to inspect all tifs
+## takes 2hrs to inspect all tifs
 library(snowfall)
 sfInit(parallel=TRUE, cpus=48)
 sfLibrary(rgdal)
@@ -481,6 +495,7 @@ sfStop()
 save.image()
 
 ## ERROR 1: /data/tt/SoilGrids250m/predicted250m/T43440/ORCDRC_M_sl1_T43440.tif, band 1: IReadBlock failed at X offset 0, Y offset 0
+#make_mosaick_ll(varn="HISTPR", in.path="/data/tt/SoilGrids250m/predicted250m", tr=cellsize, te=te, ot=metasd[which(metasd$FileName == "HISTPR_250m_ll.tif"), "DATA_FORMAT"], dstnodata=metasd[which(metasd$FileName == "HISTPR_250m_ll.tif"), "NO_DATA"], metadata=metasd[which(metasd$FileName == "HISTPR_250m_ll.tif"), sel.metasd])
 
 ## ------------- VISUALIZATIONS -----------
 
@@ -560,4 +575,3 @@ CV_NPDB_V2_ORCDRC = list(list(Observed=NPDB_V2.ovA$ORCDRC, Predicted=NPDB_V2.ovA
 save(CV_NPDB_V2_ORCDRC, file="/data/profs/SPROPS/tmp/CV_NPDB_V2_ORCDRC.rda")
 ## Plot errors:
 plot_hexbin("NPDB_V2_ORCDRC", c(0,800), main="NPDB_V2_ORCDRC", in.file="/data/profs/SPROPS/tmp/CV_NPDB_V2_ORCDRC.rda", log.plot=TRUE, out.file = "plot_CV_NPDB_V2_ORCDRC.png")
-
diff --git a/grids/models/corrupt_tif.csv b/grids/models/corrupt_tif.csv
@@ -0,0 +1 @@
+""
diff --git a/grids/models/mosaick_functions_ll.R b/grids/models/mosaick_functions_ll.R
@@ -74,13 +74,13 @@ latlon2sin = function(input.file, output.file, mod.grid, tmp.dir="/data/tmp/", p
 ## resample maps to coarser resolution ----
 aggr_SG <- function(i, r, tr=1/120, tr.metric=1000, out.dir="/data/aggregated/1km/", ti="250m", tn="1km"){
   if(missing(r)){
-    if(any(basename(i) %in% c("TAXOUSDA_250m_ll.tif", "TAXNWRB_250m_ll.tif", "TAXNWRB_300m_sin.tif", "TAXOUSDA_300m_sin.tif", paste0("TEXMHT_M_sl",1:7,"_250m_ll.tif")))){
+    if(any(basename(i) %in% c("TAXOUSDA_250m_ll.tif", "TAXNWRB_250m_ll.tif", "TAXNWRB_300m_sin.tif", "TAXOUSDA_300m_sin.tif", "GAUL_ADMIN1_landmask_300m_sin.tif", paste0("TEXMHT_M_sl",1:7,"_250m_ll.tif")))){
       r = 'near'
     } else {
       r = 'average'
     }
   }
-  if(any(basename(i) %in% c("OCSTHA_M_30cm_300m_sin.tif", "OCSTHA_M_100cm_300m_sin.tif", "OCSTHA_M_200cm_300m_sin.tif", "TAXNWRB_300m_sin.tif", "TAXOUSDA_300m_sin.tif"))){
+  if(any(basename(i) %in% c("OCSTHA_M_30cm_300m_sin.tif", "OCSTHA_M_100cm_300m_sin.tif", "OCSTHA_M_200cm_300m_sin.tif", "TAXNWRB_300m_sin.tif", "TAXOUSDA_300m_sin.tif", "GAUL_ADMIN1_landmask_300m_sin.tif"))){
     tr = tr.metric
     ti = "300m"
   }