R/convertGrid.R

In ssaxe-usgs/EASEgridR: Convert EASE-grid netCDF to WGS84

#' Convert EASE-grid netCDF to 'raster' object
#'
#' Converts an EASE-grid format netCDF file to \linkS4class{raster} object through the \code{raster} package.
#' @param nc_path Path to netCDF file. Defaults to \code{NULL}. Required if \code{lon}, \code{lat}, and \code{vals} arguments are not supplied.
#' @param name Name of variable to extract. Defaults to \code{NULL} Required if \code{lon}, \code{lat}, and \code{vals} arguments are not supplied.
#' @param lon Optional \code{vector} of longitudinal coordinates. Defaults to \code{NULL}. Required if \code{nc_path} and \code{name} arguments are not supplied.
#' @param lat Optional \code{vector} of latitudinal coordinates. Defaults to \code{NULL}. Required if \code{nc_path} and \code{name} arguments are not supplied.
#' @param vals Optional \code{matrix} of measured values. Defaults to \code{NULL}. Required if \code{nc_path} and \code{name} arguments are not supplied.
#' @param inCRS \code{CRS} object of dataset projection system. Defaults to \code{'+init=epsg:4053'}.
#' @export
#' @return \code{raster} object
#' @examples
#' convertGrid(nc_path = '.../smos.nc', name = 'Soil_Moisture')
#' convertGrid(lon = lon_vector, lat = lat_vector, vals = vals_matrix)


convertGrid <- function(nc_path = NULL,
                        name    = NULL,
                        lon     = NULL,
                        lat     = NULL,
                        vals    = NULL,
                        inCRS   = "+init=epsg:4053"){
  require(ncdf4)
  require(raster)
  require(sp)
  require(rgdal)

  # check arguments
  checkNull <- c(is.null(lon), is.null(lat), is.null(vals))
  if (all(checkNull)){ # if all lon,lat,vals args are null
    ### get values from file
    # check arguments
    if (is.null(nc_path) | is.null(name)){
      stop("If arguments 'lon', 'lat', and 'vals' are not supplied, 'nc_path' and 'name' must both be supplied.")
    }

    # open connection to netcdf
    nc <- ncdf4::nc_open(nc_path)

    # get lat, lon, and values and place into list
    sm <- list(
      x = ncdf4::ncvar_get(nc, "lon"),
      y = ncdf4::ncvar_get(nc,"lat"),
      z = ncdf4::ncvar_get(nc, name)
    )

    # close connection to netcdf
    ncdf4::nc_close(nc)

  }else if (all(!checkNull)){  # if all lon,lat,vals args are not NULL
    ### use supplied values instead of from file
    # place lat, lon, and vlaues into list
    sm <- list(
      x = lon,
      y = lat,
      z = vals
    )

  }else{
    stop("If any one of 'lon', 'lat', and 'vals' arguments are supplied, ALL must be supplied")
  }

  # create longitude SpatialPoints object
  xp <- data.frame(x = sm$x,
                   y = 0)
  sp::coordinates(xp) <- ~x+y
  sp::proj4string(xp) <- sp::CRS(inCRS)
  xp <- sp::spTransform(xp, sp::CRS("+init=epsg:3410"))

  # create latitude SpatialPoints object
  yp <- data.frame(x = 0,
                   y = sm$y)
  sp::coordinates(yp) <- ~x+y
  sp::proj4string(yp) <- sp::CRS(inCRS)
  yp <- sp::spTransform(yp, sp::CRS("+init=epsg:3410"))

  # get coordinates
  sm$xp <- sp::coordinates(xp)[,1]
  sm$yp <- sp::coordinates(yp)[,2]

  # create raster
  smr <- raster::raster(list(x = sm$xp,
                             y = sm$yp,
                             z = sm$z),
                        crs = "+init=epsg:3410")


  # return raster
  return(smr)
}