[R] Need help using GRASS within R - problem with CRS using the 'v.generalize' command

[Micha Silver]

Hello


On 08/06/2020 19:52, Loïc Valéry wrote:
> Dear all,
>
> First of all, this is my first message on the list. Therefore, please be indulgent if my message is not perfectly formatted as it should be.
>
> I am currently encountering a difficulty with GRASS 7.8 within R when using the 'v.generalize' command to smooth the contour of polygons after a segmentation step.
>
> I tried two different ways to "call" GRASS:
>
>            1 - using the RQGIS3 package

My first suggestion: no need for getting QGIS involved here. That's an 
extra layer of complication that seems unnecessary, given your goal of 
using the GRASS module, v.generalize.


>            2 - using the rgrass7 package
>
> The first method returns an error message (i.e. "proj_create_from_database: Cannot find proj.db") and therefore does not produce any result.
> The second method results in a layer of smoothed polygons but the projection reference (i.e. CRS) is lost whereas the input layer has one.
>
> Since in both cases the problem seems to be the same (i.e. GRASS fails to access the projection information), I thought it was interesting to deal with these two cases simultaneously. So, below you will find two small examples - each dealing with one of the two procedures - in order to clarify the problem.
>
> 1 - Example using the RQGIS3 package :
>
>> setwd("D:/test")
>> library(sp)
>> library(rgdal)
> rgdal: version: 1.4-8, (SVN revision 845)
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
> Path to GDAL shared files: C:/Users/toto/Documents/R/win-library/3.6/rgdal/gdal
> GDAL binary built with GEOS: TRUE
> Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
> Path to PROJ.4 shared files: C:/Users/toto/Documents/R/win-library/3.6/rgdal/proj
> Linking to sp version: 1.4-1


Please note that you have older versions of GDAL (2.2.3 here) and PROJ.4 
(4.9 here). These are currently being replaced by GDAL 3.0 and PROJ 6.3. 
You might (should) want to follow the discussion the the r-sig-geo maillist:


https://stat.ethz.ch/pipermail/r-sig-geo/2020-June/028165.html


....... (skipped all the discussion regarding RQGIS3 as I think it's not 
relevant)
>
>
> 2 - EXAMPLE USING THE RGRASS7 PACKAGE
>
> seg_poly = a SpatialPolygonsDataFrame with CRS (cf. below) :
>
>> setwd("D:/test")
>> library(rgrass7)
>>
>> # characteristics of the SpatialPolygonsDataFrame 'seg_poly' : CRS does exist
>> seg_poly
> class       : SpatialPolygonsDataFrame
> features    : 31
> extent      : 477371.3, 477397.6, 5631995, 5632020  (xmin, xmax, ymin, ymax)
> crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
> variables   : 1
> names       : Seg_ID
> min values  :      1
> max values  :     31
> Warning message:
> In proj4string(x) : CRS object has comment, which is lost in output
>> # initialization of GRASS 7.8 from R
>> initGRASS(gisBase ="C:/Program Files/GRASS GIS 7.8", home="temp/GRASS",gisDbase="temp/GRASS", use_g.dirseps.exe=F,remove_GISRC=T, override=T)
> gisdbase    temp/GRASS
> location    file19685026c56
> mapset      file196829fa7141
> rows        1
> columns     1
> north       1
> south       0
> west        0
> east        1
> nsres       1
> ewres       1
> projection  NA
>
> I suspect that the problem comes from 'projection NA' when initializing GRASS (cf. just above)


What you need to do here is setup the CRS of your new location. 
Typically, you would run initGRASS and point to a *previously created*  
LOCATION, with CRS already defined. In this case, since you are creating 
a new location, you must define it's coordinate system. (GRASS is very 
"picky" about that).

Here's a GIS Stackexchange post that explains:

https://gis.stackexchange.com/questions/183032/create-a-new-grass-database-in-r-with-crs-projection


You can derive the full proj4 string from your sp object with:


p4str = sp::proj4string(seg_poly)


Then use that to set the project parameters for the new LOCATION, with

execGRASS("g.proj", flags = "c", proj4 = p4str) Now you should be able 
to continue with...


>> execGRASS("v.generalize",flag=c("overwrite"),parameters=list(input="vec1",
> +                                                              output="GRASS_smooth_seg_poly",
> +                                                              error="GRASS_smooth_seg_poly_error",
> +                                                              method="distance_weighting",
> +                                                              threshold=1))
>
>
>
> As a novice in the field, I would be very grateful for your help.
> I remain at your entire disposal for any further information you may need to help me in finding a solution to this problem.
>
> Yours sincerely,
> Loïc
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
https://orcid.org/0000-0002-1128-1325