linkGRASS real world data usecase

Chris Reudenbach


Real world example

A typical example is the usage of an already existing project database in GRASS. GRASS organizes all data in an internal file structure that is known as gisdbase folder, a mapset and one or more locations within this mapset. All raster and vector data is stored inside this structure and the organisation is performed by GRASS. So a typical task could be to work on data sets that are already stored in an existing GRASS structure

Creating a GRASS project

Download Zensus Data

First of all we need some real world data. In this this case the gridded German 2011 micro zensus data. Download the data:

# (

of Germany. It has some nice aspects:

We also have to download the meta data description file from the above website for informations about projection and data concepts and so on.

 # we need some additional packages

# first of all we create  a project folder structure 
  link2GI::initProj(projRootDir = paste0(tempdir(),"/link2GI_examples"), 
                    projFolders =  c("run/"),
                    path_prefix = "path_",
                    global = TRUE)

# set runtime directory

# get some typical authority generated data 
 res <- curl::curl_download(url, paste0(path_run,""))

# unzip it
 unzip(res,files = grep(".csv", unzip(res,list = TRUE)$Name,value = TRUE),
       junkpaths = TRUE, overwrite = TRUE)
fn <- list.files(pattern = "[.]csv$", path = getwd(), full.names = TRUE)

Preprocessing of the data

After downloading the data we will use it for some demonstration stuff. If you have a look the data is nothing than x,y,z with assuming some projection information.

# fast read with data.table 
 xyz <- data.table::fread(paste0(path_run,"/Zensus_Bevoelkerung_100m-Gitter.csv"))


We can easy rasterize this data as it is intentionally gridded data.that means we have in at a grid size of 100 by 100 meters a value.


# clean dataframe
 xyz <- xyz[,-1]

# rasterize it according to the projection 
  r <- terra::rast(xyz, type="xyz")
 terra::crs(r) <- 3035

# map it
 p <- colorRampPalette(brewer.pal(8, "Reds"))
 # aet resolution to 1 sqkm
 mapview::mapviewOptions(mapview.maxpixels = r@ncols*r@nrows/10)
 mapview::mapview(r, col.regions = p, 
                  at = c(-1,10,25,50,100,500,1000,2500), 
                  legend = TRUE)

Setup GRASS Project

So far nothing new. Now we create a new but permanent GRASS gisbase using the spatial parameters from the raster object. As you know the linkGRASS function performs a full search for one or more than one existing GRASS installations. If a valid GRASS installation exists all parameter are setup und the package rgrass is linked.

Due to the fact that the gisdbase_exist is by default set to FALSE it will create a new structure according to the R object.

# initialize GRASS and set up a permanent structure  
link2GI::linkGRASS(x = r, 
                    gisdbase = paste0(tempdir(),"/link2GI_examples"),
                    location = "microzensus2011")   

Finally we can now import the data to the GRASS gisdbase using the rgass package functionality.

First we must convert the raster/terra object to a GeoTIFF file. Any GDAL format is possible but GeoTIFF is very common and stable.


# write it to geotiff
  terra::writeRaster(r, paste0(path_run,"/Zensus_Bevoelkerung_100m-Gitter.tif"), 
       x               overwrite = TRUE)

# import raster to GRASS

# check imported data set
                   map = "Zensus_Bevoelkerung_100m_Gitter") 

Let’s do now the same import as a vector data set. First we create a sf object. Please note this will take quite a while.

 xyz_sf = st_as_sf(xyz,
                    coords = c("x_mp_100m", "y_mp_100m"),
                    crs = 3035,
                    agr = "constant")

#map points

The GRASS gisdbase already exists. So we pass linkGRASS the argument gisdbase_exist=TRUE and import the xyz data as generic GRASS vector points.

  sf2gvec(x =  xyz_sf,
           obj_name = "Zensus_Bevoelkerung_100m_",
           gisdbase = paste0(tempdir(),"/link2GI_examples"),
           location = "microzensus2011",
           gisdbase_exist = TRUE
# check imported data set
rgrass::execGRASS('', map = "Zensus_Bevoelkerung_100m_")