inst/Kent_scripts/soilgrids_soilsgrab.R

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)
	}
bmcnellis/RSFIA documentation built on June 1, 2019, 7:40 a.m.