A typical example is the usage of an already existing project database in
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
First of all we need some real world data. In this this case the gridded German 2011 micro zensus data (see https://www.zensus2011.de/EN/Home/). 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 require(link2GI) require(curl) # 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 setwd(path_run) # get some typical authority generated data url<-"https://www.zensus2011.de/SharedDocs/Downloads/DE/Pressemitteilung/ DemografischeGrunddaten/csv_Bevoelkerung_100m_Gitter.zip; jsessionid=294313DDBB57914D6636DE373897A3F2.2_cid389?__blob=publicationFile&v=3" res <- curl::curl_download(url, paste0(path_run,"testdata.zip")) # 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)
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.
# get the filename # fast read with data.table xyz <- data.table::fread(paste0(path_run,"/Zensus_Bevoelkerung_100m-Gitter.csv")) head(xyz)
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.
require(RColorBrewer) require(raster) require(mapview) # clean dataframe xyz <- xyz[,-1] # rasterize it according to the projection r <- raster::rasterFromXYZ(xyz,crs = sp::CRS("+init=epsg: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)
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
require(link2GI) # 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
rgass7 package functionality.
First we must convert the raster object to
GeoTIFF file. Any
GDAL format is possible but
GeoTIFF is very common and stable.
require(link2GI) require(raster) require(rgrass) # write it to geotiff raster::writeRaster(r, paste0(path_run,"/Zensus_Bevoelkerung_100m-Gitter.tif"), overwrite = TRUE) # import raster to GRASS rgrass::execGRASS('r.external', flags=c('o',"overwrite","quiet"), input=paste0(path_run,"/Zensus_Bevoelkerung_100m-Gitter.tif"), output="Zensus_Bevoelkerung_100m_Gitter", band=1) # check imported data set rgrass::execGRASS('r.info', 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 sf::plot_sf(xyz_sf)
GRASS gisdbase already exists. So we pass
linkGRASS the argument
gisdbase_exist=TRUE and import the xyz data as generic GRASS vector points.
require(sf) require(sp) require(link2GI) 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('v.info', map = "Zensus_Bevoelkerung_100m_")