Thursday, 17 March 2016

r - Creating depressionless DEMs with rgrass7


I'm working with the rgrass7 package to use GRASS 7.6.0 functions within R. What I'm trying to do is the following: Import a raster DEM with a custom CRS, fill the sinks and export the depressionless raster - all via RStudio.
Some sample data can be found at: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT


I have never worked with GRASS before so I might just be missing a single command, but all my outputs are just empty rasters. Maybe this is a projection problem?


My code so far is this:


library(rgrass7)

initGRASS(gisBase = "/usr/.../grass76/",
home = tempdir(),
mapset = "PERMANENT",
override = TRUE)

# modify current mapset with custom projection
execGRASS("g.proj", flags = "c", proj4 = "+proj=aea +lat_1=25 +lat_2=50 +lat_0=37 +lon_0=87 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
execGRASS("g.mapset", flags = "c", mapset = "new_mapset")

# load DEM

execGRASS("r.in.gdal", input = "/path/to/.tif", output = "GRASS_raster", flags = c("overwrite"))

# set region (following the comment of @dmci = fixes the issue)
execGRASS("g.region", raster = raster("/path/to/.tif"))

# fill sinks
execGRASS("r.fill.dir", input = "GRASS_raster", output = "GRASS_filled_DEM", direction = "flow_test") # fill sinks
# export depressionless DEM
execGRASS("r.out.gdal", input = "GRASS_filled_DEM", output = "path/to/filled_DEM.tif")) # export filled DEM


This process does create a GRASS_raster and a filled_DEM.tif but they are just empty 1x1 raster files when I load them into RStudio.


I have tried to do this process with the GUI of GRASS. The imported DEM can be seen in the map display window, but the results of r.fill.dir are only empty rasters like this:


enter image description here



Answer



As mentioned in the comment, after importing your raster into your mapset, it is necessary to run g.region rast=youinputraster or as follows based on your script:


# set region (following the comment of @dmci = fixes the issue)
execGRASS("g.region", raster = raster("/path/to/youinputraster.tif"))

More details and options relating to g.region can be found on the GRASS documentation page


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...