I want to perform a small calculation in GRASS using the r.watershed-function. I convert my raster DEM to a SpatialGridDataFrame to import it into GRASS. However when initializing GRASS the projection information of the DEM-SpatialGridDataFrame is not kept/imported (i.e. projection(streams_g_buf) # [1] NA). Can somebody help me to spot the problem?


## get Malawi DEM
mwi_dem <- getData("alt", country="MWI", mask=TRUE)

# reproject raster to UTM 36 S
mwidem_trans <- projectRaster(mwi_dem, crs=CRS("+init=epsg:32736"))
# as SpatialGridDataFrame for GRASS
mwidem_trans2 <- as(mwidem_trans, "SpatialGridDataFrame") 

projection(mwidem_trans2) # [1] "+init=epsg:32736 +proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

## Initialize GRASS
          home=tempdir(), override=TRUE, SG=mwidem_trans2)

## Write SpatialGridDataFrame
writeRAST(mwidem_trans2, "DEM", zcol = 1, flags="overwrite")

# GRASS analyses
execGRASS("r.watershed", flags="overwrite",
          parameters=list(elevation="DEM", stream="streams_g", threshold=30))
execGRASS("r.buffer", flags="overwrite",
          parameters=list(input="streams_g", output="streams_g_buf",  distances=200))

## Read back into R
streams_g_buf <- raster(readRAST("streams_g_buf"))

projection(streams_g_buf) # [1] NA

I followed your steps up to initGRASS and then ran:

execGRASS('g.proj', flags="p")

and got:

 XY location (unprojected)

So the temp location and mapset are created without CRS information. Then when you export back R, you get NA as the CRS.

The answers here will probably help you. In short you need to create a LOCATION with proper PERMANENT mapset and CRS data, then using that LOCATION, add a mapset, pull in the R raster, then do the GRASS processing and export back to R.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.