5

How to create a new GRASS database in R with CRS projection? I am able to create a new one using the following code but it is not projected

require(raster)
require(spgrass6)
initGRASS("/Applications/GRASS-6.4.app/Contents/MacOS", gisDbase="/Users/junf/Documents",
          location='temp',mapset='PERMANENT', override=TRUE)
# gisdbase    /Users/junf/Documents 
# location    temp 
# mapset      PERMANENT 
# rows        1 
# columns     1 
# north       1 
# south       0 
# west        0 
# east        1 
# nsres       1 
# ewres       1 
# projection  NA 

Now I tried using a projected GeoTiff using the following code which is better but the projection is still NA:

r <- as(raster('/Users/junf/Documents/WatershedModel/UpperPulangi/upw_srtm30.tif'), 'SpatialGrid')
initGRASS("/Applications/GRASS-6.4.app/Contents/MacOS", gisDbase="/Users/junf/Documents", location='temp',mapset='PERMANENT', override=TRUE,SG=r)
# gisdbase    /Users/junf/Documents 
# location    temp 
# mapset      PERMANENT 
# rows        3317 
# columns     2774 
# north       950826 
# south       849428.7 
# west        690839.6 
# east        775638 
# nsres       30.56898 
# ewres       30.569 
# projection  NA 

Can someone please help me create a new GRASS database with WGS 84 (EPSG: 4326) projection or a projection that matches an existing georeferenced image for import in GRASS?

4

I came up with the following approach which is adopted from R-sig-Geo. First, create a new location via initGRASS using the default PERMANENT mapset. Next, set the desired coordinate reference system (CRS). And finally, create a new mapset which will automatically inherit the CRS from PERMANENT.

library(rgrass7)

## initialize new location and baseline mapset with required epsg code
gisBase <- system("grass72 --config path", intern = TRUE)

initGRASS(gisBase = gisBase, home = raster::tmpDir(), 
          gisDbase = "~/grassdata", location = "tmp", 
          mapset = "PERMANENT", override = TRUE)

execGRASS("g.proj", flags = "c", epsg = 4326)

## initialize new mapset inheriting projection info
execGRASS("g.mapset", flags = "c", mapset = "new_mapset")

Ideally, you should call g.mapset and g.region with print option enabled afterwards to check if everything is running smoothly.

execGRASS("g.mapset", flags = "p")
# new_mapset

execGRASS("g.region", flags = "p")
# projection: 3 (Latitude-Longitude)
# zone:       0
# datum:      wgs84
# ellipsoid:  wgs84
# north:      1N
# south:      0
# west:       0
# east:       1E
# nsres:      1
# ewres:      1
# rows:       1
# cols:       1
# cells:      1
  • Very nice and simple, this should be the best answer! Thanks. – Jun Folledo Mar 19 '17 at 8:04
3

I found that it is possible on rgrass7 but not on spgrass6: https://grasswiki.osgeo.org/wiki/GRASS_Location_Wizard

require(rgrass7)
# use epsg code
system('"C:/Program Files/GRASS GIS 7.0.3/grass70" -c epsg:32651 D:/temp/use_epsg')

# use a georeferenced raster
system('"C:/Program Files/GRASS GIS 7.0.3/grass70" -c D:/folder/srtm30.tif D:/temp/use_raster')

# connect to grass database
initGRASS(gisBase='C:/Program Files/GRASS GIS 7.0.3', gisDbase='D:/temp', location='use_raster', mapset='PERMANENT', override = TRUE)
# gisdbase    D:/temp 
# location    rtss 
# mapset      PERMANENT 
# rows        748 
# columns     944 
# north       1253702 
# south       1230844 
# west        732001.4 
# east        760855.3 
# nsres       30.55878 
# ewres       30.56557 
# projection  +proj=utm +no_defs +zone=51 +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +to_meter=1 
1

same problem stole several of my precious hours today.

It seems the only workaround is:

initGRASS(using whatever location)

Then create a new location:

execGRASS(
        'g.proj',
        parameters = list(
            georef = [your GeoTiff], 
            location = [new location name]
        )
)

initGRASS(using the new location)
1

In case you can actually use rgrass7 instead of spgrass6 you can achieve what you want with linkGRASS7 from package link2GI. I am using the development version from github here.

# devtools::install_github("gisma/link2GI")
library(link2GI)

linkGRASS7(mapview::poppendorf[[5]]) # a raster layer

gives

gisdbase    /tmp/Rtmpyd2Dh2 
location    file13e5f027f0e 
mapset      PERMANENT 
rows        254 
columns     260 
north       5515185 
south       5507565 
west        654435 
east        662235 
nsres       30 
ewres       30 
projection  +proj=utm +no_defs +zone=32 +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +to_meter=1 

Note how linkGRASS7 will by default set up your dbase and location in a temporary location (the temp folder linked to your current R session). You are able to change this by supplying additional arguments gisdbase and location.

The additional, and probably even higher benefit of link2GI is that it will automatically set the correct grass7 environment based on what it finds on your machine (therefore, linkGRAS7 is somewhat slow - you can skip that part by providing the appropriate settings directly to linkGRASS7)

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.