Loïc Valéry
2020-Jun-08 16:52 UTC
[R] Need help using GRASS within R - problem with CRS using the 'v.generalize' command
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
Micha Silver
2020-Jun-09 06:44 UTC
[R] Need help using GRASS within R - problem with CRS using the 'v.generalize' command
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 packageMy 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-1Please 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 at 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