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

Loïc Valéry |v@|ery @end|ng |rom out|ook@|r
Mon Jun 8 18:52:08 CEST 2020

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
          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 
> library(rgeos)
rgeos version: 0.5-3, (SVN revision 634)
GEOS runtime version: 3.8.0-CAPI-1.13.1 
Linking to sp version: 1.4-1 
Polygon checking: TRUE 
> library(raster)
> library(sf)
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
> library(reticulate, lib.loc="C:/Users/toto/documents/RQGIS3")
> library(RQGIS3, lib.loc="C:/Users/toto/Documents/RQGIS3")

> # set the environment to run QGIS from within R
> set_env(root="C:/Program Files/QGIS 3.12")

[1] "C:/Program Files/QGIS 3.12"

[1] "C:/Program Files/QGIS 3.12/apps/qgis"

[1] "C:/Program Files/QGIS 3.12/apps/qgis/python/plugins"

[1] "Windows"


> # Opening of the application
> open_app()
proj_create_from_database: Cannot find proj.db
proj_create_from_database: Cannot find proj.db

However, this does not completely prevent the package from working because I can find the function I want, call it and enter the input parameters. Here is an excerpt :

> # Search for the desired function within QGIS 3.12
> find_algorithms(search_term="generalize", name_only=TRUE)
[1] "grass7:v.generalize"

> # Information on how to use the selected function (here, the function 'generalize' from GRASS 7 within QGIS 3.12)
> get_usage(alg="grass7:v.generalize")
v.generalize (grass7:v.generalize)
Vector based generalization.
Input parameters

input: Input layer
         Parameter type:  QgsProcessingParameterFeatureSource
         Accepted data types:
                 - str: layer ID
                 - str: layer name
                 - str: layer source
                 - QgsProcessingFeatureSourceDefinition
                 - QgsProperty
                 - QgsVectorLayer
type: Input feature type
         Parameter type:  QgsProcessingParameterEnum
         Available values:
                 - 0: line
                 - 1: boundary
                 - 2: area
         Accepted data types:
                 - int
                 - str: as string representation of int
e.g. 1
                 - QgsProperty

> # Indicates the mandatory parameters
> params <- get_args_man(alg = "grass7:v.generalize")
Choosing default values for following parameters:
type: 0
method: 0
See get_options('grass7:v.generalize') for all available options.

> # Input of mandatory parameters
> params$input<-"D:/test/seg_poly.shp"
> params$type<-"0"
> params$method<-"7"
> params$threshold<-1
> params$output<-"D:/test/smooth_seg_poly.shp"
> params$error<- "D:/test/smooth_seg_poly_error.shp"

But when I execute the function, R returns an error message that is related to the previous error message (i.e. when I executed the 'open_app()' command). Please find below the second error message :

> # Execution of the selected function
> out <- run_qgis(alg = "grass7:v.generalize",
+                 params = params,
+                 load_output = TRUE)

ERROR 1: PROJ: proj_identify: Cannot find proj.db
proj_create_from_wkt: Cannot find proj.db
proj_identify: Cannot find proj.db
Error in py_call_impl(callable, dots$args, dots$keywords) : 

  QgsProcessingException: There were errors executing the algorithm.

Detailed traceback: 
  File "C:/Program Files/QGIS 3.12/apps/qgis/python/plugins\processing\tools\general.py", line 106, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "C:/Program Files/QGIS 3.12/apps/qgis/python/plugins\processing\core\Processing.py", line 181, in runAlgorithm
    raise QgsProcessingException(msg)


To help you find the solution to the problem, here are the results of some commands specifying the working environment :

> Sys.info()
       sysname        release        version       nodename        machine          login           user effective_user 
     "Windows"      "8.1 x64"   "build 9600"     "BERNACHE"       "x86-64"         "toto"         "toto"         "toto" 

> sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)
Matrix products: default

[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252    LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    
attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] RQGIS3_1.0.1.9000 reticulate_1.16   sf_0.9-3          raster_3.1-5      rgeos_0.5-3       rgdal_1.4-8       sp_1.4-2         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6       compiler_3.6.3     pillar_1.4.4       class_7.3-15       tools_3.6.3        jsonlite_1.6.1     tibble_3.0.1      
 [8] lifecycle_0.2.0    lattice_0.20-38    pkgconfig_2.0.3    rlang_0.4.6        Matrix_1.2-18      DBI_1.1.0          cli_2.0.2         
[15] rstudioapi_0.11    parallel_3.6.3     e1071_1.7-3        stringr_1.4.0      vctrs_0.3.0        hms_0.5.3          classInt_0.4-3    
[22] grid_3.6.3         glue_1.4.1         R6_2.4.1           fansi_0.4.1        readr_1.3.1        magrittr_1.5       codetools_0.2-16  
[29] ellipsis_0.3.1     units_0.6-6        assertthat_0.2.1   KernSmooth_2.23-16 stringi_1.4.6      crayon_1.3.4      

> qgis_session_info()

[1] "3.0.4"


[1] "3.12.0-București"

[1] "2.3.2"

 Finally, since the problem probably comes from accessing the proj.db file, I ran the command Sys.getenv(“PROJ_LIB”) and R finds two locations…

> Sys.getenv("PROJ_LIB")
[1] "C:/Program Files/QGIS 3.12/share/proj;C:/Users/toto/Documents/R/win-library/3.6/rgdal/proj"

Please, note that, when searching on my computer, the file proj.db is also found in the following locations (which R does not seem to find) :
1.       C:/Users/toto/Documents/R/win-library/3.6/sf/proj
2.       C:/Programmes/GRASS GIS 7.8/share/proj
3.       C:/Programmes/QGIS 3.12/apps/Python3.7/lib/site-packages/pyproj/proj-dir/share/proj


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)

Running GRASS : everything seems to be going fine...

> # specifies to GRASS that objects of type "sp" are used
> rgrass7::use_sp()
> # export of the vector in GRASS format from the sp object 'seg_poly'
> writeVECT(seg_poly,"vec1",v.in.ogr_flags=c("o", "overwrite"), driver="ESRI Shapefile")

> # run the 'v.generalize' command to smooth the contours of polygons
> 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))

...but when retrieving the result, the CRS has been lost :

> # retrieving the result of the 'v.generalize' command in a SpatialPolygonDataFrame : 'smooth_seg_poly'
> smooth_seg_poly<-readVECT("GRASS_smooth_seg_poly", with_prj=T, driver="ESRI Shapefile")
Exporting 31 areas (may take some time)...
VOUTOG~1 terminé. 31 features (Polygon type) written to <GRASS_sm>
(ESRI_Shapefile format).
OGR data source with driver: ESRI Shapefile 
Source: "D:\test\temp\GRASS\file19685026c56\file196829fa7141\.tmp\unknown", layer: "GRASS_sm"
with 31 features
It has 2 fields
Warning message:
In vColumns(vname) : vColumns: v.info -c output not in two columns:
Displaying column types/names for database connection of layer <1>:
> # characteristics of the SpatialPolygonsDataFrame 'smoooth_seg_poly' : CRS no longer exists !!!
> smooth_seg_poly
class       : SpatialPolygonsDataFrame 
features    : 31 
extent      : 477371.3, 477397.6, 5631995, 5632020  (xmin, xmax, ymin, ymax)
crs         : NA 
variables   : 2
names       : cat, Seg_ID 
min values  :   1,      1 
max values  :  31,     31 

If it helps some info on my R session and on some environment variables (i.e. PROJ_LIB and GRASS_PROJSHARE)

> sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 8.1 x64 (build 9600)

Matrix products: default

[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252    LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    

attached base packages:
[1] tools     stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] RSAGA_1.3.0         plyr_1.8.6          shapefiles_0.7      foreign_0.8-75      gstat_2.0-6         exactextractr_0.4.0 sf_0.9-3           
 [8] rgrass7_0.2-1       RSQLite_2.2.0       XML_3.99-0.3        RQGIS3_1.0.1.9000   reticulate_1.16     smoothr_0.1.2       maptools_1.0-1     
[15] link2GI_0.4.3       glue_1.4.1          listviewer_3.0.0    raster_3.1-5        rgeos_0.5-3         rgdal_1.5-8         sp_1.4-2           

loaded via a namespace (and not attached):
 [1] pkgload_1.0.2      bit64_0.9-7        jsonlite_1.6.1     R.utils_2.9.2      assertthat_0.2.1   blob_1.2.1         remotes_2.1.1     
 [8] sessioninfo_1.1.1  pillar_1.4.4       backports_1.1.7    lattice_0.20-38    digest_0.6.25      R.oo_1.23.0        htmltools_0.4.0   
[15] Matrix_1.2-18      pkgconfig_2.0.3    devtools_2.3.0     purrr_0.3.4        intervals_0.15.2   processx_3.4.2     tibble_3.0.1      
[22] usethis_1.6.1      ellipsis_0.3.1     withr_2.2.0        cli_2.0.2          magrittr_1.5       crayon_1.3.4       memoise_1.1.0     
[29] ps_1.3.3           R.methodsS3_1.8.0  fs_1.4.1           fansi_0.4.1        xts_0.12-0         xml2_1.3.2         class_7.3-15      
[36] FNN_1.1.3          pkgbuild_1.0.8     prettyunits_1.1.1  hms_0.5.3          lifecycle_0.2.0    stringr_1.4.0      callr_3.4.3       
[43] compiler_3.6.3     e1071_1.7-3        spacetime_1.2-3    rlang_0.4.6        classInt_0.4-3     units_0.6-6        grid_3.6.3        
[50] rstudioapi_0.11    htmlwidgets_1.5.1  testthat_2.3.2     codetools_0.2-16   DBI_1.1.0          roxygen2_7.1.0     R6_2.4.1          
[57] zoo_1.8-8          knitr_1.28         bit_1.1-15.2       rprojroot_1.3-2    KernSmooth_2.23-16 readr_1.3.1        desc_1.2.0        
[64] stringi_1.4.6      parallel_3.6.3     Rcpp_1.0.4.6       vctrs_0.3.0        xfun_0.14         
> Sys.getenv("PROJ_LIB")
[1] "C:/Users/toto/Documents/R/win-library/3.6/rgdal/proj"
> Sys.getenv("GRASS_PROJSHARE")
[1] "NA\\share\\proj"

I guess that the "NA\\share\\proj" just above is not normal but I don't know precisely what path should be the right one...  and this anomaly may not be the only problem...

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,

More information about the R-help mailing list