2

I'm trying to calculate subglacial sink volume by using the rgrass7 package to utilize GRASS functions in R. However, the sinks found by the function seem to be very small: enter image description here

When I run the similar functions in ArcGIS, the result makes much more sense as there are all kinds of sink sizes: enter image description here

(These two images are only examples, I have verified the same results with multiple glacier rasters.)

Are there noteworthy differences between the two routines? Does GRASS not fill sinks to the pour point? Can I somehow give more parameters to r.fill.dir like the fill height etc.?

My code is:

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

# modify current mapset with 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")

### run GRASS commands 
execGRASS("r.in.gdal", input = myraster, output = paste0("GRASS_raster"),
          flags = c("overwrite"))
execGRASS("g.region", raster = "GRASS_raster")
execGRASS("r.fill.dir", input = "GRASS_raster", output = "GRASS_filled_DEM", direction = "flow_test", flags = c("overwrite")) # fill sinks
execGRASS("r.out.gdal", input = "GRASS_filled_DEM", output = paste0(output, "filled_DEM.tif"), flags = c("overwrite"))
1

I found out a possible answer (not sure, if it could be applied universally, but for me it worked). r.fill.dir needs to be executed multiple times in order to fill every unresolved area to the pour point.

### run GRASS commands 
execGRASS("r.in.gdal", input = myraster, output = paste0("GRASS_raster"),
          flags = c("overwrite"))
execGRASS("g.region", raster = "GRASS_raster")
execGRASS("r.fill.dir", input = "GRASS_raster", output = "GRASS_filled_DEM", direction = "flow_test", flags = c("overwrite")) # fill sinks
execGRASS("r.fill.dir", input = "GRASS_filled_DEM", output = "GRASS_filled_DEM", direction = "flow_test", flags = c("overwrite")) # repeat until every area is resolved
execGRASS("r.out.gdal", input = "GRASS_filled_DEM", output = paste0(output, "filled_DEM.tif"), flags = c("overwrite"))
| improve this answer | |

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.