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

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Loïc Valéry
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")

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

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

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

$platform
[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
GRASS_OUTPUT_TYPE_PARAMETER: 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"
> params$GRASS_OUTPUT_TYPE_PARAMETER<-"0"

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

locale:
[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()

$gdal
[1] "3.0.4"

$grass7
[1] FALSE

$qgis_version
[1] "3.12.0-București"

$saga
[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



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)


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)...
 100%
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

locale:
[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,
Loïc
______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

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

Micha Silver-2
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
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.