setwd("C:/Work/CABBI/R")
#sorgdata <- read.table("C:/Work/CABBI/R/esorg_obsdata.tsv",header=T,sep="\t")
sorgdata <- read.table("C:/Work/CABBI/R/esorg_validdata.tsv",header=T,sep="\t")
sorgdata$longlat <- paste(sorgdata$latitude,sorgdata$longitude,sep="_")
wthfiles <- sorgdata[!duplicated(sorgdata$longlat),]
wthfiles <- wthfiles[!is.na(wthfiles$longitude),]
urls <- paste0("https://rest.soilgrids.org/query?lon=",wthfiles$longitude,"&lat=",wthfiles$latitude)
## query urls one at a time, skipping NA centroids
for (i in 1:length(urls)) #loop through the vector of urls (aka queries)
{
download.file(url=urls[i],
destfile=paste0("soils/",wthfiles[i,"SITE"],"_",wthfiles[i,"longlat"],".json")) #download.file is an R utility that sends the query and directs the response to a file on disk, giving it the supplied filename
}
#############################################################################
#### ONCE YOU HAVE THE JSON FILES, below code extracts and formats the data into soils.in files
#### SOILS.IN FIELDS:
# Column 1 - Minimum depth of soil layer (cm)
# Column 2 - Maximum depth of soil layer (cm)
#Column 3 - Bulk density of soil layer (g/cm^3) === BLDFIE gives BD in kg/m^3
# Column 4 - Field capacity of soil layer, volumetric === AWCh1-3 give volumetric water capacity (fraction)
# Column 5 - Wilting point of soil layer, volumetric === WWP gives AWC until wilting point (vol fraction)
## Column 6 - Evaporation coefficient for soil layer (currently not being used)
## Column 7 - Percentage of roots in soil layer, these values must sum to 1.0
# Column 8 - Fraction of sand in soil layer, 0.0 - 1.0 === SNDPPT gives sand content 1-100
# Column 9 - Fraction of clay in soil layer, 0.0 - 1.0 === CLYPPT gives clay content 1-100
# Column 10 - Organic matter in soil layer, fraction 0.0 - 1.0 === ORCDRC gives SOC in g/kg of soil, 7 layers
## Column 11 - Minimum volumetric soil water content below wilting point for soil
# layer, soil water content will not be allowed to drop below this
# value
# Column 12 - Saturated hydraulic conductivity of soil layer in centimeters per
# second === fn "ksat" calculates it according to the Saxton equation in units of mm/h
# Column 13 - pH of soil layer === PHIHOX gives pH of soil in water
## JS SOIL LAYER POINTDEPTHS: 0, 5, 15, 30, 60, 100, 200
#### READ EACH VARIABLE INTO ITS OWN NAMED VECTOR AND RECAST VALUES FROM POINT MEASUREMENTS TO DEPTH-INTERVAL AVERAGES
setwd("C:/Work/CABBI/R")
library(jsonlite)
source("idemaz_rcropmod.R")
myjss <- list.files("soils",pattern="json")
jslayers <- c(0, 5, 15, 30, 60, 100, 200) #specific depths at which measurements were taken
jsmindep <- round(jslayers-(c(0,diff(jslayers)/2)),0) #center the soils.in layers on the depths at which json point-measurements were taken
jsmaxdep <- c(jsmindep[2:7],250)
spacepad <- function(x)
{
pad <- strrep(" ",times=10-nchar(x))
return(paste0(x,pad))
}
for (i in 1:length(myjss))
{
myjs <- read_json(paste0("soils/",myjss[i])) #read current file into an R list
thickn1 <- jsmaxdep - jsmindep
mindep2 <- jsmindep
maxdep3 <- jsmaxdep
blkden4 <- unname(unlist(myjs$properties$BLDFIE[[2]]))/1000 #given as kg/m3, convert to g/cm3
fldcap5 <- unname(unlist(myjs$properties$AWCh1[[2]]))/100 + unname(unlist(myjs$properties$WWP[[2]]))/100 #given as 0-100 percent, convert to 0-1 fraction
wltpnt6 <- unname(unlist(myjs$properties$WWP[[2]]))/100 #""
evpcoe7 <- c(0.8,0.2,rep(0,5)) #not in use, pad with missing value flags
pctrts8 <- c(0.2,0.3,0.3,0.2,0,0,0) #guesstimate for now
sndfrc9 <- unname(unlist(myjs$properties$SNDPPT[[2]]))/100 #""
clyfrc10 <- unname(unlist(myjs$properties$CLYPPT[[2]]))/100 #""
orgmt11 <- unname(unlist(myjs$properties$ORCDRC[[2]]))/1000 #given in g/kg, convert to g/g fraction 0-1
mnswc12 <- c(0.08,0.08,0.05,0.02,0,0,0) ## not a standard soil property, for now filling with values from a stock soils.in file
stcnd13 <- ksat(soc=pmin(orgmt10*100,8),sand=sndfrc8,clay=clyfrc9)/(3600*10) ## saturated hydr conductivity, cm/s, use idemaz/Rcropr implementation of Saxton 2006 equations; convert result from mm/hr to cm/s
stcnd13 <- max(0.0001,signif(stcnd12,5)) #ksat returns a big mess of digits so round it
phh2o14 <- unname(unlist(myjs$properties$PHIHOX[[2]]))/10 # given as ph*10 convert to pH
## cbind the columns and write the files
soilsin <- cbind(thickn1,mindep2,maxdep3,blkden4,fldcap5,wltpnt6,evpcoe7,pctrts8, #bind the vectors as columns into a matrix
sndfrc9,clyfrc10,orgmt11,mnswc12,stcnd13,phh2o14)
## apply to each column
soilsin <- rbind(soilsin[1,],soilsin[1,],soilsin[2:7,]) #split first layer, keep ID values (DayCent wants 0-2, 2-5)
soilsin[2,1:3] <- c(3,2,5) #put in the new layer thickness/min/maxdeps
soilsin[3,1:2] <- c(5,5) #""
soilsin <- rbind(soilsin[1:5,],soilsin[5:6,],soilsin[6:7,],soilsin[7,],soilsin[7,],soilsin[7:8,])
soilsin[,"mindep2"] <- c(0,2,5,10,20,30,45,60,75,90,105,120,150)
soilsin[,"maxdep3"] <- c(2,5,10,20,30,45,60,75,90,105,120,150,180)
soilsinchar <- apply(soilsin,c(1,2),formatC,format="f",digits=4)
soilsinchar <- soilsinchar[,2:ncol(soilsinchar)]
write.table(soilsinchar,file=paste0("soils/",sub("json","soils.in",myjss[i])),sep="\t",col.names=F,row.names=F,quote=F)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.