rm(list=ls(all=T)) gc() table.hwsd <- read.csv("HWSD_DATA.csv") # add some extension lines to describe our self-made continental sea types lastentry <- length(table.hwsd[,1]) lastmugl <- max(as.integer(table.hwsd[,"MU_GLOBAL"])) if (lastmugl >= 32100) { print("ERROR: last MU_GLOBAL in HWSD_DATA.csv >= 32100 already"); quit() } lastmugl <- 32099 # first add possible values to the SU_SYM74 column # here we insert our self-made type symbols into column SU_SYM74, they are then propagated to column SU_SYM90 below levels(table.hwsd[,"SU_SYM74"]) <- c(levels(table.hwsd[,"SU_SYM74"]), "CSGL", "CSAS", "CSRS", "CSGA", "CSPG", "CSEC", "CSDS", "CSGI") # Great Lakes lastentry <- lastentry+1 mugl_gl <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSGL", 0) # Aral Sea lastentry <- lastentry+1 mugl_as <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSAS", 0) # Red Sea lastentry <- lastentry+1 mugl_rs <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSRS", 0) # Gulf of Aden lastentry <- lastentry+1 mugl_ga <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSGA", 0) # Persian Gulf lastentry <- lastentry+1 mugl_pg <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSPG", 0) # English Channel lastentry <- lastentry+1 mugl_ec <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSEC", 0) # Denmark Street lastentry <- lastentry+1 mugl_ds <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSDS", 0) # Gibraltar lastentry <- lastentry+1 mugl_gi <- lastmugl <- lastmugl+1 table.hwsd[lastentry,c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_SYM74","ADD_PROP")] <- c(lastentry, lastmugl, 0, 100.00, "CSGI", 0) # replace missing values with 0 table.hwsd[which(is.na(table.hwsd[,"T_GRAVEL"])),"T_GRAVEL"] <- 0 table.hwsd[which(is.na(table.hwsd[,"S_GRAVEL"])),"S_GRAVEL"] <- 0 table.hwsd[which(table.hwsd[,"SU_SYM90"]=="HSf" & is.na(table.hwsd[,"T_SAND"])),c("T_SAND","T_SILT","T_CLAY")] <- 0 # fill missing values in column SU_SYM90 with values from column SU_SYM74 # note that this also concerns greenland glaciers with MU_GLOBAL=6998,SU_SYM74="GG",SU_SUM90="" levels(table.hwsd[,"SU_SYM90"]) <- c(levels(table.hwsd[,"SU_SYM90"]), levels(table.hwsd[,"SU_SYM74"])) table.hwsd[which(is.na(table.hwsd[,"SU_CODE90"])),"SU_SYM90"] <- table.hwsd[which(is.na(table.hwsd[,"SU_CODE90"])),"SU_SYM74"] table.hwsd[,"SHARE"] <- as.double(table.hwsd[,"SHARE"])/100 # length(table.hwsd[which(is.na(table.hwsd[,"T_SILT"])),"SU_SYM90"]) # # summary(table.hwsd[which(table.hwsd[,"SU_SYM90"]=="DS" | table.hwsd[,"SU_SYM90"]=="FP" | table.hwsd[,"SU_SYM90"]=="GG" | table.hwsd[,"SU_SYM90"]=="HD" | table.hwsd[,"SU_SYM90"]=="IS" | table.hwsd[,"SU_SYM90"]=="NI" | table.hwsd[,"SU_SYM90"]=="RK" | table.hwsd[,"SU_SYM90"]=="ST" | table.hwsd[,"SU_SYM90"]=="UR" | table.hwsd[,"SU_SYM90"]=="WR"),c("T_SAND","T_SILT","T_CLAY","S_SAND","S_SILT","S_CLAY")]) # # unique(table.hwsd[which(table.hwsd[,"T_SAND"]==0 & table.hwsd[,"T_SILT"]==0 & table.hwsd[,"T_CLAY"]==0),c("SU_SYM74","SU_SYM85","SU_SYM90")]) imu.null <- which(table.hwsd[,"T_SAND"]==0 & table.hwsd[,"T_SILT"]==0 & table.hwsd[,"T_CLAY"]==0) imu.na <- which(is.na(table.hwsd[,"T_SAND"])) file.hwsd <- file("hwsd.bil","rb") frac.soil <- array(data=0,dim=c(720,360)) frac.soil_ts <- array(data=0,dim=c(720,360)) frac.sea <- array(data=0,dim=c(720,360)) frac.null <- array(data=0,dim=c(720,360)) frac.999 <- array(data=0,dim=c(720,360)) frac.DS <- array(data=0,dim=c(720,360)) # DS Dunes & shifting sands frac.FP <- array(data=0,dim=c(720,360)) # FP Fishpond frac.GG <- array(data=0,dim=c(720,360)) # GG Glaciers & permanent snow frac.HD <- array(data=0,dim=c(720,360)) # HD Humanly disturbed frac.IS <- array(data=0,dim=c(720,360)) # IS Island frac.NI <- array(data=0,dim=c(720,360)) # NI No data frac.RK <- array(data=0,dim=c(720,360)) # RK Rock debris frac.ST <- array(data=0,dim=c(720,360)) # ST Salt flats frac.UR <- array(data=0,dim=c(720,360)) # UR Urban frac.WR <- array(data=0,dim=c(720,360)) # WR Inland water frac.CSGL <- array(data=0,dim=c(720,360)) # CSGL Continental Sea: Great Lakes #frac.CSAS <- array(data=0,dim=c(720,360)) # CSAS Continental Sea: Aral Sea #frac.CSRS <- array(data=0,dim=c(720,360)) # CSRS Continental Sea: Red Sea #frac.CSGA <- array(data=0,dim=c(720,360)) # CSGA Continental Sea: Gulf of Aden #frac.CSPG <- array(data=0,dim=c(720,360)) # CSPG Continental Sea: Persian Gulf #frac.CSEC <- array(data=0,dim=c(720,360)) # CSEC Continental Sea: English Channel #frac.CSDS <- array(data=0,dim=c(720,360)) # CSDS Continental Sea: Denmark Street #frac.CSGI <- array(data=0,dim=c(720,360)) # CSGI Continental Sea: Gibraltar #summary(frac.soil+frac.sea+frac.null+frac.999+frac.DS+frac.FP+frac.GG+frac.HD+frac.IS+frac.NI+frac.RK+frac.ST+frac.UR+frac.WR) data.t_gravel <- array(data=0,dim=c(720,360)) data.t_sand <- array(data=0,dim=c(720,360)) data.t_silt <- array(data=0,dim=c(720,360)) data.t_clay <- array(data=0,dim=c(720,360)) data.t_oc <- array(data=0,dim=c(720,360)) data.t_gravel_ts <- array(data=0,dim=c(720,360)) data.t_sand_ts <- array(data=0,dim=c(720,360)) data.t_silt_ts <- array(data=0,dim=c(720,360)) data.t_clay_ts <- array(data=0,dim=c(720,360)) data.t_oc_ts <- array(data=0,dim=c(720,360)) data.s_gravel <- array(data=0,dim=c(720,360)) data.s_sand <- array(data=0,dim=c(720,360)) data.s_silt <- array(data=0,dim=c(720,360)) data.s_clay <- array(data=0,dim=c(720,360)) data.s_oc <- array(data=0,dim=c(720,360)) chunk.hwsd <- array(dim=c(60,43200)) # a function to replace 0-entries in chunk.hwsd inside a given bounding box with a given # value mu_global. the lat and lon values are indices in the global hwsd map of 43200x21600 # points. ilat is the number of the currently loaded latitude band of 60 rows. patch_hwsd <- function(lonmin, lonmax, latmin, latmax, ilat, mu_global) { #print(paste(lonmin, lonmax, latmin, latmax, ilat, mu_global)) latmin <- latmin - (ilat-1)*60 latmax <- latmax - (ilat-1)*60 if (latmax < 1) return() # entire rectangle above band if (latmin > 60) return() # entire rectangle below band if (latmin < 1) latmin <- 1 if (latmax > 60) latmax <- 60 print(paste(latmin, latmax)) chunk.hwsd[which(chunk.hwsd[latmin:latmax,lonmin:lonmax] == 0)] <<- mu_global } for(ilat in 1:360) { print(paste("Process latitude band ",ilat,sep="")) seek(file.hwsd,where=(ilat-1)*60*43200*2,origin="start") for(i in 1:60) chunk.hwsd[i,] <- readBin(file.hwsd,integer(),size=2,n=43200) # patch 0-values with above self-defined values for special regions # Great Lakes patch_hwsd(10540, 12780, 4900, 5850, ilat, mugl_gl) patch_hwsd(12780, 14500, 4750, 5360, ilat, mugl_gl) # Aral Sea patch_hwsd(28580, 29020, 5170, 5600, ilat, mugl_as) # Red Sea patch_hwsd(25480, 26830, 7100, 9300, ilat, mugl_rs) # Gulf of Aden patch_hwsd(26700, 27750, 8970, 9560, ilat, mugl_ga) # Persian Gulf patch_hwsd(27335, 27890, 7140, 7925, ilat, mugl_pg) patch_hwsd(27890, 28420, 7476, 7925, ilat, mugl_pg) patch_hwsd(28420, 28640, 7540, 7980, ilat, mugl_pg) patch_hwsd(28640, 28775, 7720, 8100, ilat, mugl_pg) # English Channel patch_hwsd(20990, 22000, 4540, 5000, ilat, mugl_ec) # Denmark Street patch_hwsd(22560, 23400, 3600, 4380, ilat, mugl_ds) # Gibraltar patch_hwsd(20700, 21100, 6350, 6666, ilat, mugl_gi) for(ilon in 1:720) { cell <- chunk.hwsd[,(1:60)+(ilon-1)*60] mu <- unique(cell[1:3600]) for(n in 1:length(mu)) { mufrac <- length(which(cell==mu[n]))/3600 if(mu[n]==0) frac.sea[ilon,ilat] <- mufrac else if(mu[n]==999) frac.999[ilon,ilat] <- mufrac else { imu <- which(table.hwsd[,"MU_GLOBAL"]==mu[n]) for(i in 1:length(imu)) { if(is.element(imu[i],imu.null)) frac.null[ilon,ilat] <- frac.null[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"] else if(is.element(imu[i],imu.na)) switch(as.character(table.hwsd[imu[i],"SU_SYM90"]), "DS" = frac.DS[ilon,ilat] <- frac.DS[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "FP" = frac.FP[ilon,ilat] <- frac.FP[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "GG" = frac.GG[ilon,ilat] <- frac.GG[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "HD" = frac.HD[ilon,ilat] <- frac.HD[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "IS" = frac.IS[ilon,ilat] <- frac.IS[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "NI" = frac.NI[ilon,ilat] <- frac.NI[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "RK" = frac.RK[ilon,ilat] <- frac.RK[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "ST" = frac.ST[ilon,ilat] <- frac.ST[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "UR" = frac.UR[ilon,ilat] <- frac.UR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], # for now, treat all continental sea values just as water "CSGL" = frac.CSGL[ilon,ilat] <- frac.CSGL[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSAS" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSRS" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSGA" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSPG" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSEC" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSDS" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "CSGI" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"], "WR" = frac.WR[ilon,ilat] <- frac.WR[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]) else { frac.soil[ilon,ilat] <- frac.soil[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"] data.t_gravel[ilon,ilat] <- data.t_gravel[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_GRAVEL"] data.t_sand[ilon,ilat] <- data.t_sand[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_SAND"] data.t_silt[ilon,ilat] <- data.t_silt[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_SILT"] data.t_clay[ilon,ilat] <- data.t_clay[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_CLAY"] data.t_oc[ilon,ilat] <- data.t_oc[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_OC"] if(is.finite(table.hwsd[imu[i],"S_SAND"])) { frac.soil_ts[ilon,ilat] <- frac.soil_ts[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"] data.t_gravel_ts[ilon,ilat] <- data.t_gravel_ts[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_GRAVEL"] data.t_sand_ts[ilon,ilat] <- data.t_sand_ts[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_SAND"] data.t_silt_ts[ilon,ilat] <- data.t_silt_ts[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_SILT"] data.t_clay_ts[ilon,ilat] <- data.t_clay_ts[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_CLAY"] data.t_oc_ts[ilon,ilat] <- data.t_oc_ts[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"T_OC"] data.s_gravel[ilon,ilat] <- data.s_gravel[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"S_GRAVEL"] data.s_sand[ilon,ilat] <- data.s_sand[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"S_SAND"] data.s_silt[ilon,ilat] <- data.s_silt[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"S_SILT"] data.s_clay[ilon,ilat] <- data.s_clay[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"S_CLAY"] data.s_oc[ilon,ilat] <- data.s_oc[ilon,ilat] + mufrac*table.hwsd[imu[i],"SHARE"]*table.hwsd[imu[i],"S_OC"] } } } } } } # print(frac.GG[,ilat]) # print(frac.GG[270:280,40:45]) } #print(frac.GG[275,42]) # 60 0.01*1 # 1972 0.09805556*1 # 1443 0.27*(0.5+0.3+0.1+0.1) # 1544 0.05916667*1 # 1572 0.4408333*(0.6+0.3+0.1) # 1276 0.1213889*(0.5+0.5) # 1455 0.0005555556*(0.9+0.1) # # 0.01*1+0.09805556*1+0.27*(0.5+0.3+0.1+0.1)+0.05916667*1+0.4408333*(0.6+0.3+0.1)+0.1213889*(0.5+0.5)+0.0005555556*(0.9+0.1) close(file.hwsd) #summary(frac.soil[,13]+frac.sea[,13]+frac.null[,13]+frac.DS[,13]+frac.FP[,13]+frac.GG[,13]+frac.HD[,13]+frac.IS[,13]+frac.NI[,13]+frac.RK[,13]+frac.ST[,13]+frac.UR[,13]+frac.WR[,13]) #summary(frac.soil+frac.sea+frac.null+frac.999+frac.DS+frac.FP+frac.GG+frac.HD+frac.IS+frac.NI+frac.RK+frac.ST+frac.UR+frac.WR) #summary(data.t_sand+data.t_silt+data.t_clay) #summary(data.s_sand+data.s_silt+data.s_clay) save.image("process_hwsd1.RData") frac.soil <- frac.soil[,360:1] frac.soil_ts <- frac.soil_ts[,360:1] frac.sea <- frac.sea[,360:1] frac.null <- frac.null[,360:1] frac.DS <- frac.DS[,360:1] # DS Dunes & shifting sands frac.FP <- frac.FP[,360:1] # FP Fishpond frac.GG <- frac.GG[,360:1] # GG Glaciers & permanent snow frac.HD <- frac.HD[,360:1] # HD Humanly disturbed frac.IS <- frac.IS[,360:1] # IS Island frac.NI <- frac.NI[,360:1] # NI No data frac.RK <- frac.RK[,360:1] # RK Rock debris frac.ST <- frac.ST[,360:1] # ST Salt flats frac.UR <- frac.UR[,360:1] # UR Urban frac.CSGL <- frac.CSGL[,360:1] # CSGL Continental Sea: Great Lakes frac.WR <- frac.WR[,360:1] # WR Inland water #print(frac.GG[270:280,40:45]) #print(frac.GG[270:280,(361-40):(361-45)]) data.t_gravel <- data.t_gravel[,360:1]/frac.soil data.t_sand <- data.t_sand[,360:1]/frac.soil data.t_silt <- data.t_silt[,360:1]/frac.soil data.t_clay <- data.t_clay[,360:1]/frac.soil data.t_oc <- data.t_oc[,360:1]/frac.soil data.t_gravel_ts <- data.t_gravel_ts[,360:1]/frac.soil_ts data.t_sand_ts <- data.t_sand_ts[,360:1]/frac.soil_ts data.t_silt_ts <- data.t_silt_ts[,360:1]/frac.soil_ts data.t_clay_ts <- data.t_clay_ts[,360:1]/frac.soil_ts data.t_oc_ts <- data.t_oc_ts[,360:1]/frac.soil_ts data.s_gravel <- data.s_gravel[,360:1]/frac.soil_ts data.s_sand <- data.s_sand[,360:1]/frac.soil_ts data.s_silt <- data.s_silt[,360:1]/frac.soil_ts data.s_clay <- data.s_clay[,360:1]/frac.soil_ts data.s_oc <- data.s_oc[,360:1]/frac.soil_ts save.image("process_hwsd2.RData") # # data.hwsd <- data.hwsd[,360:1] # # library(fields) # image.plot(x=seq(-179.75,179.75,len=720), y=seq(-89.75,89.75,len=360), data.hwsd, zlim=c(1,max(data.hwsd)), xlim=c(-180,180), ylim=c(-90,90), col=rainbow(n=50,start=0,end=5/6), xlab="",ylab="")