R/magpie2raster.R

In pik-piam/geodata: Landuse geodata library

### script to convert magpie object to raster

## check magpie2raster(ratio[,,"crop"]) has two data dimensions but is not recogniced so
## CHECK SHOULD LOOK AT ALL DATA DIMENSIONS!!

## to do:
# enable an option that several magpie dimensions can also become a stack not a list of rasters
# 




#' Convert MAgPIE object to raster
#' 
#' MAgPIE object with a resolution of 0.5x0.5 deg and 59199 cells is converted
#' to raster data. Data with several years will be converted to a RasterStack,
#' data with several data dimensions to a list of RasterLayers or RasterStacks.
#' 
#' 
#' @usage magpie2raster(x)
#' @param x A MAgPIE object
#' @return A RasterLayer/-Stack or a list of such
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{
#' a <- new.magpie(cells_and_regions = paste(magclass:::magclassdata$half_deg$region,1:59199,sep="."),
#'                 years = c(2000,2005), names = c("crop", "past"), fill = rnorm(59199*4))
#' r <- magpie2raster(a)
#' r
#' }
#' 
#' @export magpie2raster
#' @importFrom magclass is.magpie fulldim getYears getNames
#' @importFrom raster res<- stack res
magpie2raster <- function(x) {
  
  if(!is.magpie(x)) stop("x has to be a magpie object")
  if(fulldim(x)[[1]][1]!=59199) stop("x has to be on 0.5 deg resolution, and contain 59199 cells")
  
  

  # the index of the magpie cells in the raster
  rast <- raster::raster()
  raster::res(rast) <- 0.5
  coord <- as.matrix(ludata::getCoordinates(degree = T))
  index <- raster::cellFromXY(rast, coord)
  
  # easiest case: magpie object has only one year and dimension, then is is just created as RasterLayer
  if ( all(fulldim(x)[[1]][-1]==1) ) {
    rast <- raster::raster()
    raster::res(rast) <- 0.5
    rast[index] <- as.vector(x)
    names(rast) <- getYears(x, as.integer = F)
  
  # magpie object has several years, then is is  created as RasterStack
  } else if ( fulldim(x)[[1]][2] > 1 & all(fulldim(x)[[1]][c(-1,-2)] == 1) ) {
    rast <-  raster::stack()
    
    for (y in getYears(x, as.integer = F)) {
      ylayer <- raster::raster()
      res(ylayer) <- 0.5
      ylayer[index] <- as.vector(x[,y,])
      names(ylayer) <- y
      rast <- stack(rast, ylayer)
    }
  
  # magpie object has several data dimensions but only one year, then a list of RasterLayers is created
  } else if ( fulldim(x)[[1]][2] == 1 &  any(fulldim(x)[[1]][c(-1,-2)] > 1) ) {
    rast <- list()
    
    for (n in getNames(x)) {
      nlayer <- raster::raster()
      raster::res(nlayer) <- 0.5
      nlayer[index] <- as.vector(x[,,n])
      names(nlayer) <- getYears(x, as.integer = F)
      rast[n] <- nlayer
    }
    
  # magpie object has several data dimensions and several year dimensions, then a list of RasterStacks is created
  } else if ( fulldim(x)[[1]][2] > 1 & any(fulldim(x)[[1]][c(-1,-2)] > 1) ) {
    rast <- list()
    
    for (n in getNames(x)) {
      
      yrast <-  raster::stack()
      for (y in getYears(x, as.integer = F)) {
        ylayer <- raster::raster()
        res(ylayer) <- 0.5
        ylayer[index] <- as.vector(x[,y,n])
        names(ylayer) <- y
        yrast <- stack(yrast, ylayer)
      }
      rast[n] <- yrast
    }
    
  }
    
  return(rast)  
  
}