RSAGA geoprocessing instructions - ncss-tech/geo-pit GitHub Wiki

This document displays some R batch functions for generating DEM derivatives using the RSAGA R package. It is assumes the reader has already compiled a series of DEM following the nedInstrucitons document.

SAGA is an open-source GIS that was originally developed in 1996 as a terrain analysis toolbox, known as DiGem. Since then it has become a fully fledged GIS, with additional modules for vector geoprocessing, image analysis, and geostatistics. While not as well documented at GRASS or ArcGIS it offers an intuitive interface, and includes a range algorithms not found elsewhere. Through the use of the RSAGA package, SAGA modules can be called from R, and workflows developed. Unlike other GIS, SAGA utilizes significant RAM instead of using file caching. This makes SAGA fast, but it can also overwhelm a computer if to many large rasters are loaded. However I've noticed when using a solid state drive (SSD) I can process rasters than exceded my 16GB of RAM for certain SAGA modules that only use small local neighorhoods.

To begin, the necessary libaries must be loaded, as well as the custom batch functions.

library(gdalUtils)
library(RSAGA)

source("C:/Users/Stephen/Documents/Github/geoprocessing/gdalUtilsFunctions.R")
source("C:/Users/Stephen/Documents/Github/geoprocessing/rsagaFunctions.R")

Next the proper GDAL and RSAGA path has to be set. The first GDAL location is the default path on my work computer, the second my personal computer. If this isn't set gdalUtils will do a brute force search of your computer, which usually finds GDAL 1.7 instead of the GDAL 10.1. The new version has additional features, which many these batch functions use.

gdal_setInstallation(search_path = "C:/ProgramData/QGIS/QGISDufour/bin", rescan = TRUE)
gdal_setInstallation(search_path = "C:/OSGeo4W64/bin", rescan = T)
myenv <- rsaga.env(path="C:/Program Files (x86)/SAGA-GIS")

Next numerous parameters need to be set which get used later by many of the functions or commands. Modify these file paths and lists as necessary. For example, I organized my files by "C:/geodata/project_data/11ATL"", so 11 will have to replace by 10 or 2 for your respective Regions.

office <- c("8VIC", "8PHO")
sdat.p <- paste0("M:/geodata/project_data/", office, "/sdat/")

# Generater raster names
res10 <- paste0(sdat.p, "ned10m_", office)
res30 <- paste0(sdat.p, "ned30m_", office)
radiusD <- 2
radiusV <- round(1000/30/2-1, 0)

g10 <- list(
  slopeR = paste0(res10, "_slopeR", 1+2*radiusD),
  slope  = paste0(res10, "_slope", 1+2*radiusD),
  slopeD = paste0(res10, "_slopeD", 1+2*radiusD),
  aspect = paste0(res10, "_aspect", 1+2*radiusD),
  cupro  = paste0(res10, "_cupro", 1+2*radiusD),
  cucon  = paste0(res10, "_cucon", 1+2*radiusD),
  cutan  = paste0(res10, "_cutan", 1+2*radiusD),
  cumax  = paste0(res10, "_cumax", 1+2*radiusD),
  cumin  = paste0(res10, "_cumin", 1+2*radiusD)
  )
g30 <- list(
  elev     = paste0(res30, "_elev", 1+2*radiusD),
  slope    = paste0(res30, "_slope", 1+2*radiusD),
  slopeR   = paste0(res30, "_slopeR", 1+2*radiusD),
  slopeD   = paste0(res30, "_sloped", 1+2*radiusD),
  aspect   = paste0(res30, "_aspect", 1+2*radiusD),
  mvalleys  = paste0(res30, "_mvalleys"),
  msummits  = paste0(res30, "_msummits"),
  caarea   = paste0(res30, "_caarea"),
  caheight = paste0(res30, "_cheight"),
  wetness  = paste0(res30, "_wetness"),
  strahler = paste0(res30, "_strahler"),
  z2stream = paste0(res30, "_z2stream")
  )

Convert GTiff to SAGA

dem10.tif  <- lapply(strsplit(res10, "/sdat"), paste0, collapse="")
dem10.tif  <- paste0(dem10.tif, ".tif")
dem30.tif  <- lapply(strsplit(res30, "/sdat"), paste0, collapse="")
dem30.tif  <- paste0(dem30.tif, ".tif")
dem10.sdat <- paste0(res10, ".sdat")
dem30.sdat <- paste0(res30, ".sdat")

gdal_GTiff2SAGA(dem10.tif, dem10.sdat)
gdal_GTiff2SAGA(dem30.tif, dem30.sdat)

Calculate local derivatives

dem10 <- paste0(res10, ".sgrd")
attach(lapply(g10, paste0, ".sgrd"))

rsaga.d1(dem10, 2)
rsaga.d2(dem10, 2)
rsaga.d3(dem10, 2)
# Convert radians to percent, degrees=radians*180/pi, 180/pi=57.29578 
rsaga.grid.calculus(slopeD, slopeD, slope, "tan(a*(1/57.29578))*100")
# Rescale curvatures so they can be exported as UInt16 to save file size
rsaga.grid.calculus(cupro, cupro, cupro, "10000*a")
rsaga.grid.calculus(cucon, cucon, cucon, "-10000*a")
rsaga.grid.calculus(cumin, cumin, cumin, "10000*a")
rsaga.grid.calculus(cumax, cumax, cumax, "10000*a")
rsaga.grid.calculus(cucon, slopeD, cutan, "a*sin(b/57.29578)")

dem30 <- paste0(res30, ".sgrd")
attach(lapply(g30, paste0, ".sgrd"))

rsaga.d0(dem30, 2)
rsaga.d1(dem30, 2)
rsaga.grid.calculus(slopeD, slopeD, slope, "tan(a*(1/57.29578))*100")
rsaga.grid.calculus(slopeD, slopeD, slopeR, "a*(1/57.29578)")
rsaga.mrvbf(dem30, mvalleys, msummits)

Create a hydrologically correct DEM

# Create a regional mosaic
elev.sdat <- paste0(g30$elev, ".sdat")
mosaicList(list(elev.sdat), "E:/geodata/project_data/REGION11/ned30m_R11_elev5.tif")


# Create a copy of elev5 to mask---------------------------------------------------------
gdal_translate(
  src_dataset = "E:/geodata/project_data/REGION11/ned30m_R11_elev5.tif",
  dst_dataset = "E:/geodata/project_data/REGION11/ned30m_R11_elev5_masked.tif",
  overwrite = TRUE,
  verbose = TRUE
)


# Extract the water bodies into shapefiles and reproject-----------------------------------
state <- c("NV", "CA")
nhd    <- paste0("D:/geodata/hydrography/NHDH_", state, ".gdb")
nhd_wb <- paste0("D:/geodata/hydrography/NHDH_", state, "_wb.shp")

for(i in seq(nhd)){
  ogr2ogr(
    src_datasource_name = nhd[i],
    dst_datasource_name = nhd_wb[i],
    layer = "NHDWaterbody",
    t_srs = "EPSG:5070",
    overwrite = TRUE,
    verbose = TRUE,
    progress = TRUE)
}


# Mask the water bodies from the Regional DEM------------------------------------------------
# Seems to take exceptionally long for the States touching the Great Lakes. If done separately in OSGeo4W Shell or SAGA you can monitor their progress.
for(i in seq(nhd_wb)){
  cat(paste(format(Sys.time(), "%Y-%m-%d %H:%M:%S"), "burning", nhd_wb[i], "\n"))
  gdal_rasterize(
    src_datasource = paste0("E:/geodata/hydrography"),
    dst_filename = "E:/geodata/project_data/REGION11/ned30m_R11_elev5_masked.tif",
    l = paste0("NHDH_", state[i], "_wb"),
    where = "AreaSqKm > 0.04047",
    b = 1,
    burn = -99999,
    verbose = TRUE
  ) 
}


# Subset the masked Regional DEM into MLRA office subsets----------------------------------
dem30.sdat<- paste0(res30, "_elev5_masked.sgrd")
batchSubsetSAGA("E:/geodata/project_data/REGION11/ned30m_R11_elev5_masked.tif", dem30.sdat, nlcdpath)


# Calculate hydrologcal derivatives (wetness index and relief above streams)---------------
dem <- paste0(g30$elev, "_masked.sgrd")

rsaga.fill(dem)

dem <- paste0(g30$elev, "_masked_filled.sgrd")

rsaga.ca(dem)
rsaga.twi(slopeR, caarea)
rsaga.strahler(dem, strahler, 4)
rsaga.reclassify(strahler, -3, -1)

streams <- paste0(g30$elev, "_strahler_rc0.sgrd")

rsaga.ofd(dem, streams)


# Write SAGA to GTiff--------------------------------------------------------------------------
int16 <- c("slope", "aspect")
in16.sdat <- lapply(g10[int16], paste0, ".sdat")
int16.tif <- unlist(lapply(strsplit(int16.sdat, "/sdat"), paste0, collapse=""))
int16.tif <- paste0(strsplit(int16.tif, ".sdat"), ".tif")
gdal_SAGA2GTiff(int16.sdat, int16.tif, "Int16", -32768)

slopeshape <- paste0("E:/geodata/project_data/ned10m_", office, "_slopeshape.tif")
int16 <- c("cupro", "cutan", "slope")
int16.sdat <- lapply(g10[int16], paste0, ".sdat")
gdal.stack(int16.sdat, slopeshape, "Int16", -32768)

flt <- c("z2stream", "mvalleys", "wetness")
flt.sdat <- lapply(g30[flt], paste0, ".sdat")
flt.tif <- unlist(lapply(strsplit(flt.sdat, "/sdat"), paste0, collapse=""))
flt.tif <- paste0(strsplit(flt.tif, ".sdat"), ".tif")
gdal_SAGA2GTiff(flt.sdat, flt.tif, "Float32", -99999)

mosaicList(g10$slope, "E:/geodata/project_data/11REGION/ned10m_11R_slope5.tif", "Int16", c("COMPRESS=DEFLATE", "TILED=YES", "BIGTIFF=YES"), -32768)