[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")
$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
More information about the R-help
mailing list