Chapter 11 SNODAS

11.1 How to use SNODAS

SNODAS data are available from going to the web page, clicking ‘Go to FTP,’ and downloading the files of your choosing. The masked and unmasked data files are organized on the FTP site in separate directories labeled masked and unmasked. Within these two directories are subdirectories labeled by a 4-digit year. Within the year directories, there are subdirectories for the months of the year of the form MM_mon where MM is the two-digit month number and mon is the three-character month abbreviation. Each month directory contains the tarred archive file, usually one for each day of the month. Each day has each variable separated into different files. The user guide has a table displaying the variable to product code conversions.

Below, we illustrate how to access SNODAS data using the NCAR/rwrfhydro package in R.

11.2 SNODAS Code example

Retrieving SWE and Snow Depth

library(rwrfhydro)

# getPoints ----------------------------------------------------------------
# Edited from rwrfhydro/GetSnodasPointTs, which had outdated variable names that caused an error
getPoints <- function(bDatePOSIXct,eDatePOSIXct,snodasDir,lat,lon,quiet=TRUE){
  #First, calculate SNODAS lat/lon coordinates
  snodasCoords <- CalcSnodasCoords()
  latMin <- snodasCoords$Lat[which.min(snodasCoords$Lat)]
  lonMin <- snodasCoords$Lon[which.min(snodasCoords$Lon)]
  latMax <- snodasCoords$Lat[which.max(snodasCoords$Lat)]
  lonMax <- snodasCoords$Lon[which.max(snodasCoords$Lon)]
  res <- 0.00833333333333300
  
  #Sanity check on lat/lon
  if((lat < latMin) | (lat > latMax)){
    warning("Error: Provided latitude should range from 0-90.")
    return(0)
  }
  if((lon < lonMin) | (lon > lonMax)){
    warning("Error: Provided longitude should range from -180.0 to 0.0")
    return(0)
  }
  
  #Perform date analysis to determine time difference between beginning ending dates
  dUnits <- "days"
  diff1 <- difftime(eDatePOSIXct,bDatePOSIXct,units = dUnits)
  nSteps <- diff1 <- as.numeric(diff1)
  dSec <- diff1*24*3600
  dt <- 24*3600
  
  data <- data.frame()
  
  #Calculate SNODAS x,y coordinates using lat/lon pair
  rowInd <- floor((lat - (latMin - (res/2.0)))/res) + 1
  colInd <- floor((lon - (lonMin - (res/2.0)))/res) + 1
  
  for (step in 0:(nSteps)){
    #establish date of current time step
    dCurrent <- bDatePOSIXct + dt*step
    
    if(quiet == FALSE){
      print(paste0('Extracting SNODAS for: ',strftime(dCurrent,"%Y-%m-%d")))
    }
    
    #Establish #SNODAS file name
    snodasFile <- paste0(snodasDir,'/SNODAS_',strftime(dCurrent,"%Y%m%d"),'.nc')
    
    if(!file.exists(snodasFile)){
      warning('Error: SNODAS file: ',snodasFile,' not found.')
      return(0)
    }
    
    #Open NetCDF file
    nc <- ncdf4::nc_open(snodasFile)
    #Sanity check
    varNames <- names(nc$var)
    if((varNames[1] != "SWE") & (varNames[2] != "snowDepth")){
      warning("Error: Unexpected variables found in: ",snodasFile)
      return(0)
    }
    
    #Extract fill value
    fillValueSWE <- ncdf4::ncatt_get(nc, varid="SWE", attname="_FillValue")
    fillValueSD <- ncdf4::ncatt_get(nc, varid="snowDepth", attname="_FillValue")
    
    if(fillValueSWE[1] != TRUE){
      warning("Error: FillValue for SWE not found.")
      return(0)
    }
    if(fillValueSD[1] != TRUE){
      warning("Error: FillValue for snow depth not found.")
      return(0)
    }
    #Pull data out
    sweData <- ncdf4::ncvar_get(nc, varid="SWE", start=c(colInd,rowInd,1), count=c(1,1,1))
    sdData <- ncdf4::ncvar_get(nc, varid="snowDepth", start=c(colInd,rowInd,1), count=c(1,1,1))
    
    #Check for missing values
    if(sweData == fillValueSWE[2]){
      sweData <- NA
    }else{
      #Convert to mm
      sweData <- sweData
    }
    if(sdData == fillValueSD[2]){
      sdData <- NA
    }else{
      #Convert to mm
      sdData <- sdData
    }
    
    #Close NetCDF file
    ncdf4::nc_close(nc)
    
    #Assign values to a temporary data frame
    #Note SNEQV is snow water equivalent and SNOWH is snow depth. This is to 
    #match the LSM output format. 
    dfTemp <- data.frame(POSIXct = dCurrent, SNEQV = sweData, SNOWH = sdData, 
                         units = "mm")
    
    #Merge with existing data frame
    data <- plyr::rbind.fill(data,dfTemp)
  } # end for time loop
  return(data)
} 

# Pull SNODAS for each day used in analysis
for (i in 1:31){
  if(i < 10) date = paste0(0,i)
  else date = i
  
  # January
  snodasGot <- GetSnodasDepthSweDate(as.POSIXct(paste0('2017-01-',date)),
                                outputDir = path.expand('/Users/ameyer/RShiny_Microclim/Data/SNODAS'))
  if(snodasGot) snodasList <- ReadSnodasDepthSweDate(as.POSIXct(paste0('2017-01-',date)),
                                outputDir = path.expand('/Users/ameyer/RShiny_Microclim/Data/SNODAS'))
  PutSnodasNcdf(snodasList)

  # July
  snodasGot <- GetSnodasDepthSweDate(as.POSIXct(paste0('2017-07-',date)),
                                     outputDir = path.expand('/Users/ameyer/RShiny_Microclim/Data/SNODAS'))
  if(snodasGot) snodasList <- ReadSnodasDepthSweDate(as.POSIXct(paste0('2017-07-',date)),
                                                     outputDir = path.expand('/Users/ameyer/RShiny_Microclim/Data/SNODAS'))
  PutSnodasNcdf(snodasList, outputDir = path.expand('/Users/ameyer/RShiny_Microclim/Data/SNODAS'))
  
  # Clean up files
  unlink(path.expand('/Users/ameyer/RShiny_Microclim/Data/SNODAS/*.dat.gz'))
}


# Get dataframes for each location for both months
snodasJanuaryWA <- getPoints(as.POSIXct('2017-01-01'),as.POSIXct('2017-01-31'),
                             snodasDir=getwd(),lat=47,lon=-118.57)
snodasJulyWA <- getPoints(as.POSIXct('2017-07-01'),as.POSIXct('2017-07-31'),
                             snodasDir=getwd(),lat=47,lon=-118.57)
snodasJanuaryCO <- getPoints(as.POSIXct('2017-01-01'),as.POSIXct('2017-01-31'),
                             snodasDir=getwd(),lat=40.87,lon=-104.73)
snodasJulyCO <- getPoints(as.POSIXct('2017-07-01'),as.POSIXct('2017-07-31'),
                             snodasDir=getwd(),lat=40.87,lon=-104.73)