From ac5cc85c670e0fbdc5e6f4f63bdc6f24f8316437 Mon Sep 17 00:00:00 2001 From: rmartinezferia Date: Fri, 23 Jun 2017 14:13:36 -0500 Subject: [PATCH] Added download example for RD --- ApsimSoil-SSURGO.Rproj | 13 ----- _code/downloadDataBase.R | 106 --------------------------------------- _code/downloadExample.R | 18 +++++++ _code/download_ssurgo.R | 46 ++++++----------- 4 files changed, 34 insertions(+), 149 deletions(-) delete mode 100644 ApsimSoil-SSURGO.Rproj delete mode 100644 _code/downloadDataBase.R create mode 100644 _code/downloadExample.R diff --git a/ApsimSoil-SSURGO.Rproj b/ApsimSoil-SSURGO.Rproj deleted file mode 100644 index 8e3c2eb..0000000 --- a/ApsimSoil-SSURGO.Rproj +++ /dev/null @@ -1,13 +0,0 @@ -Version: 1.0 - -RestoreWorkspace: Default -SaveWorkspace: Default -AlwaysSaveHistory: Default - -EnableCodeIndexing: Yes -UseSpacesForTab: Yes -NumSpacesForTab: 2 -Encoding: UTF-8 - -RnwWeave: Sweave -LaTeX: pdfLaTeX diff --git a/_code/downloadDataBase.R b/_code/downloadDataBase.R deleted file mode 100644 index 1398457..0000000 --- a/_code/downloadDataBase.R +++ /dev/null @@ -1,106 +0,0 @@ -#install.packages("FedData") -library(FedData) -source("_code/SaxtonRawls.R") - -# Inputs ########################################### - -SiteName <- "SorensonFarm" -FieldExtent <- c(-93.740682, -93.741712, 42.01000, 42.010759) - -# Download SSURGO data from extent ################# - -a <- polygon_from_extent(raster::extent(FieldExtent[1],FieldExtent[2],FieldExtent[3],FieldExtent[4]), - proj4string="+proj=longlat") - -x <- get_ssurgo(template=a, label=SiteName) - -# Extract useful data ############################# - -component <- x$tabular$component -chorizon <- x$tabular$chorizon -mapunit <-x$tabular$mapunit -majcomp <- component[component$majcompflag == "Yes", c("compname", - "taxclname", - "drainagecl", - "mukey","cokey","slope.r","hydgrp")] - -### Plot mapunits -#plot(x$spatial) -#invisible(text(getSpPPolygonsLabptSlots(x$spatial), labels=as.character((majcomp$compname)), cex=1)) - -# Merge data ######################################## - -majcomp<- merge(data.frame(mukey=as.character(x$spatial$MUKEY), - musym=as.character(x$spatial$MUSYM)), - majcomp) - - -# Calculate % of area -area <- data.frame(musym=x$spatial$MUSYM, area=NA) -for (i in 1:length(x$spatial$MUSYM)) area$area[i] <- x$spatial@polygons[[i]]@area -area$area <- area$area/sum(area$area) - -# Merge into dataset -majcomp <- merge(majcomp, area) -majcomp <- merge(majcomp,chorizon) - -# Create horizons -h <- majcomp[order(majcomp$compname, majcomp$hzdept.r),] -h <- h[,c("compname","area","slope.r", "hydgrp", - "hzdept.r","hzdepb.r","hzthk.r", - "sandtotal.r","claytotal.r", - "dbthirdbar.r","dbovendry.r", - "om.r","ksat.r", - "wfifteenbar.r","wthirdbar.r","ph1to1h2o.r")] - -names(h) <- c("name","area","slope","hydrogroup","top","bottom","thick","sand","clay","wetbd","drybd", "om","ksat","ll","dul","ph") - -# Create new variables ########################### - -## Center -h$center <- h$top + h$thick/2 - -# Runoff - -hydgrp <- data.frame(matrix(data=c(61,64,68,71,73,76,80,83,81,84,88,91,84,87,91,94), - ncol = 4,nrow = 4, - dimnames = list(c("0-2","2-5","5-10","10-100"), - LETTERS[1:4]))) - - -hydgrp$"A/B" <- 0.5*(hydgrp[,1]+hydgrp[,1+1]) -hydgrp$"A/C" <- 0.5*(hydgrp[,1]+hydgrp[,1+2]) -hydgrp$"A/D" <- 0.5*(hydgrp[,1]+hydgrp[,1+3]) - -hydgrp$"B/C" <- 0.5*(hydgrp[,2]+hydgrp[,2+1]) -hydgrp$"B/D" <- 0.5*(hydgrp[,2]+hydgrp[,2+2]) - -hydgrp$"C/D" <- 0.5*(hydgrp[,3]+hydgrp[,3+1]) - -require(tidyr) -hydgrp$key <-row.names(hydgrp) - -hydgrp %>% - group_by(key) %>% - gather(key="hydgrp",value="CN2",-key) - - -for(i in length(names(hydgrp))){ -# while(i < length(names(hydgrp)) ) - - 0.5*(hydgrp[,1]+hydgrp[,1+1]) - -} - - - -if your soil has a slope of 2% and belongs to Hydrological group of B, then the CN2 value is 73 - 5 = 68 -if your soil has a slope of 2% and belongs to Hydrological group of B/D, then you calculate 2 CN2 values and make the average, here is 73-5=68 and 84-5=79, so (68+79)/2 = 73.5 - - -#h$bd <- SaxtonRawls(pSand=h$sand ,pClay=h$clay, pOM=h$om)$BD -h$bd <- ifelse( h$wetbd < 0.9, 0.9, ifelse(h$wetbd > 1.8, 1.8,h$wetbd)) - -#h$ll.sr <- SaxtonRawls(pSand=h$sand ,pClay=h$clay, pOM=h$om)$LL15 -#h$dul.sr <- SaxtonRawls(pSand=h$sand ,pClay=h$clay, pOM=h$om)$DUL -h$sat <- SaxtonRawls(pSand=h$sand ,pClay=h$clay, pOM=h$om)$SAT diff --git a/_code/downloadExample.R b/_code/downloadExample.R new file mode 100644 index 0000000..1c02914 --- /dev/null +++ b/_code/downloadExample.R @@ -0,0 +1,18 @@ +# Example of function that downloads the SSURGO database for a field extent + +#install.packages("FedData") +require(tidyverse) +source("_code/download_ssurgo.R") + +# (down)load ssurgo data +h <- download_SSURGO(SiteName = "Sorenson", + # Set soil layer structure + soilLayer_breaks = c(5,20,50,80,120,180), + # Set field extent (in Decimal Degrees) + north=42.010759, # Latitude + south=42.01000, # Latitude + east=-93.740682, # Longitude + west=-93.741712) # Longitude + +# Access the downloaded data +h$soils \ No newline at end of file diff --git a/_code/download_ssurgo.R b/_code/download_ssurgo.R index 5f722f4..53ac6c2 100644 --- a/_code/download_ssurgo.R +++ b/_code/download_ssurgo.R @@ -14,9 +14,6 @@ download_SSURGO <- function(SiteName = "field2", soilLayer_breaks = c(1,5,10,20,30,50,80,120,160,200), interpolation_Method = c("linear", "constant")) { - - #SiteName = "field1";north=42.010759;south=42.01000;east=-93.740682;west=-93.741712;map=T - # Load packages require(sp) require(FedData) @@ -27,7 +24,7 @@ download_SSURGO <- function(SiteName = "field2", proj4string="+proj=longlat") x <- get_ssurgo(template=a, label=SiteName) - + # Extract useful data ############################# @@ -105,13 +102,13 @@ download_SSURGO <- function(SiteName = "field2", #soil$ll.sr <- SaxtonRawls(pSand=soil$sand ,pClay=soil$clay, pOM=soil$om)$LL15 #soil$dul.sr <- SaxtonRawls(pSand=soil$sand ,pClay=soil$clay, pOM=soil$om)$DUL - soil$ksat <- soil$ksat*1000/1.157 # mm/day + soil$ksat <- soil$ksat*100/1.157 # mm/day soil$sat <- SaxtonRawls(pSand=soil$sand ,pClay=soil$clay, pOM=soil$om)$SAT/100 soil$ph <- 0.52+1.06*soil$ph #pH 1:5 - soil$oc <- soil$om/1.72 # % + soil$OC <- soil$om/1.72 # % soil$U <- ifelse(soil$clay<=20,5+0.175*soil$clay, ifelse(soil$clay<=40,7.5+0.05*soil$clay, @@ -119,21 +116,22 @@ download_SSURGO <- function(SiteName = "field2", ifelse(soil$clay<=70,12.75-0.075*soil$clay, ifelse(soil$clay<=80,11-0.05*soil$clay,0))))) # mm - soil$conna <- ifelse(soil$clay<=30,0.025*soil$clay+3.25, + soil$cona <- ifelse(soil$clay<=30,0.025*soil$clay+3.25, ifelse(soil$clay<=50,4, ifelse(soil$clay<=70,-0.025*soil$clay+5.25, ifelse(soil$clay<=80,3.5,0)))) # mm/d^5 soil$PO <- 1-soil$bd/2.65 - soil$salb <- 0.13 - soil$FINERT <- ifelse(soil$center<=1,0.4, + soil$Salb <- 0.13 + + soil$FInert <- ifelse(soil$center<=1,0.4, ifelse(soil$center<=10,0.4, ifelse(soil$center<60,0.008*soil$center+0.32, ifelse(soil$center<=120,0.8, ifelse(soil$center<180,0.0032*soil$center+0.42, ifelse(soil$center<=300,0.99,0)))))) #(0-1) - soil$FBIOM <- ifelse(soil$center<=10,0.04, + soil$FBiom <- ifelse(soil$center<=10,0.04, ifelse(soil$center<=20,0.055-0.0015*soil$center, ifelse(soil$center<=30,0.03-0.0005*soil$center, ifelse(soil$center<60,0.0216-0.0002*soil$center, @@ -149,9 +147,9 @@ download_SSURGO <- function(SiteName = "field2", soil$DiffusSlope <- 16 - soil$CnRed <- 20 #residue runoff + soil$CNRed <- 20 #residue runoff - soil$CnCov <- 0.8 + soil$CNCov <- 0.8 soil$EnrAcoeff <- 7.4 @@ -171,29 +169,16 @@ download_SSURGO <- function(SiteName = "field2", soil$AirDry <- ifelse(soil$center<=15,0.5,ifelse(soil$center<=30,0.9,1))*soil$ll + # Bin by layer soil$layer <- .bincode(soil$center, breaks=c(-1,soilLayer_breaks[soilLayer_breaks>0])) soil %>% + select(-hydrogroup) %>% group_by(layer) %>% - summarise(top=min(center), - bottom=max(center), - area=mean(area), - sand=mean(sand), - clay=mean(clay), - bd=mean(bd), - ll=mean(ll), - dul=mean(dul), - ksat=mean(ksat), - sat=mean(sat), - oc=mean(oc), - ph=mean(ph), - FBIOM=mean(FBIOM), - FINERT=mean(FINERT), - U = mean(U) - conna = mean(conna) - ) -> soil + mutate(thick = (max(center)-min(center))*10) %>% + summarise_all("mean") -> soil # Save df into list ########### @@ -205,7 +190,8 @@ download_SSURGO <- function(SiteName = "field2", names(soils) <- (unique(h$name)) - return(soils) + return(list(coordinates = c(mean(north,south),mean(east,west)), + soils = soils)) } ### Example