gmeta() returns a proj4string

Should getLocationProj() (used by gmeta() and maybe others) switch to WKT? GRASS can return the latter with the g.proj -w command.

Thanks, perhaps, with an argument to restore legacy behaviour. When was the -w flag introduced?

When was the -w flag introduced?

No idea. I din't find it at https://trac.osgeo.org/grass.

We need to know which version introduced it in order to keep work reproducible looking backwards. It is not well documented, seems to have been in 7.8 and 7.6. The website upgrade has made things prettier but much harder to use. Probably 7.6 is a reasonable cutoff, but please check.

When was the -w flag introduced?

No idea. I din't find it at https://trac.osgeo.org/grass.

GRASS development has moved to github one year ago or so: https://github.com/OSGeo/grass. The full history of each command can be explored there or accessing through the direct link at the bottom of every manual page. For g.proj:

the history is here: https://github.com/OSGeo/grass/commits/master/general/g.proj

I've seen some recent commits on WTK by @metzm, but I'm not sure if that's what needed here

Thank you @veroandreo .

The flag appears in the documentation at least since 27 Apr 2008 (7.0), it may perhaps go further back. See documentation file https://github.com/OSGeo/grass/blob/db49180dd772ddb7643b633d57ad00c58914a519/general/g.proj/description.html.

This is the oldest commit of this github repo.

Example from here:

Convert the PROJ.4 projection description contained in a text file to WKT format:

cat proj4.description | g.proj -w proj4=-

Unfortunately, this doesn't give guidance, because -wef has existed for a long time (probably GRASS 6) to emit the contents of an ESRI Shapefile .prj file. I guess 7.6 is a reasonable cutoff for WKT2_2019, which is the dialect that is interesting here.

Please fork and install from source - the projection component will be WKT for GRASS >= 7.6 unless the user sets an optional argument to FALSE.

Thank you for this feature Roger, I will try it.

WKT will have a central role in GRASS 8's CRS handling it seems, supporting WKT1 & WKT2, see https://github.com/OSGeo/grass/pull/976 and more pull requests in https://github.com/OSGeo/grass/milestone/4?closed=1.

Ignore WKT1, use only WKT2_2019 if possible, like PROJ. Report back when you have tried it, that is ASAP.

It works perfectly here. :+1:

library(rgrass7)
#> Loading required package: XML
#> GRASS GIS interface loaded with GRASS version: (GRASS not running)
rgdal::rgdal_extSoftVersion()
#>           GDAL GDAL_with_GEOS           PROJ             sp 
#>        "3.1.3"         "TRUE"        "7.1.1"        "1.4-4"
initGRASS(gisBase = link2GI::findGRASS()[1, 1],
          mapset = "PERMANENT",
          home = tempdir())
#> gisdbase    /tmp/RtmpL4Y03f 
#> location    file5c2a32bba3d7 
#> mapset      PERMANENT 
#> rows        1 
#> columns     1 
#> north       1 
#> south       0 
#> west        0 
#> east        1 
#> nsres       1 
#> ewres       1 
#> projection:
#>  XY location (unprojected)
execGRASS("g.proj", "c", epsg = 31370)
#> Default region was updated to the new projection, but if you have multiple
#> mapsets `g.region -d` should be run in each to update the region from the
#> default
#> Projection information updated
execGRASS("g.proj", "w")
#> PROJCS["Belge 1972 / Belgian Lambert 72",
#>     GEOGCS["Belge 1972",
#>         DATUM["Reseau_National_Belge_1972",
#>             SPHEROID["International 1924",6378388,297,
#>                 AUTHORITY["EPSG","7022"]],
#>             TOWGS84[-99.059,53.322,-112.486,0.419,-0.83,1.885,-1],
#>             AUTHORITY["EPSG","6313"]],
#>         PRIMEM["Greenwich",0,
#>             AUTHORITY["EPSG","8901"]],
#>         UNIT["degree",0.0174532925199433,
#>             AUTHORITY["EPSG","9122"]],
#>         AUTHORITY["EPSG","4313"]],
#>     PROJECTION["Lambert_Conformal_Conic_2SP"],
#>     PARAMETER["latitude_of_origin",90],
#>     PARAMETER["central_meridian",4.36748666666667],
#>     PARAMETER["standard_parallel_1",51.1666672333333],
#>     PARAMETER["standard_parallel_2",49.8333339],
#>     PARAMETER["false_easting",150000.013],
#>     PARAMETER["false_northing",5400088.438],
#>     UNIT["metre",1,
#>         AUTHORITY["EPSG","9001"]],
#>     AXIS["Easting",EAST],
#>     AXIS["Northing",NORTH],
#>     AUTHORITY["EPSG","31370"]]
gmeta()
#> gisdbase    /tmp/RtmpL4Y03f 
#> location    file5c2a32bba3d7 
#> mapset      PERMANENT 
#> rows        1 
#> columns     1 
#> north       1 
#> south       0 
#> west        0 
#> east        1 
#> nsres       1 
#> ewres       1 
#> projection:
#>  PROJCS["Belge 1972 / Belgian Lambert 72",
#>     GEOGCS["Belge 1972",
#>         DATUM["Reseau_National_Belge_1972",
#>             SPHEROID["International 1924",6378388,297,
#>                 AUTHORITY["EPSG","7022"]],
#>             TOWGS84[-99.059,53.322,-112.486,0.419,-0.83,1.885,-1],
#>             AUTHORITY["EPSG","6313"]],
#>         PRIMEM["Greenwich",0,
#>             AUTHORITY["EPSG","8901"]],
#>         UNIT["degree",0.0174532925199433,
#>             AUTHORITY["EPSG","9122"]],
#>         AUTHORITY["EPSG","4313"]],
#>     PROJECTION["Lambert_Conformal_Conic_2SP"],
#>     PARAMETER["latitude_of_origin",90],
#>     PARAMETER["central_meridian",4.36748666666667],
#>     PARAMETER["standard_parallel_1",51.1666672333333],
#>     PARAMETER["standard_parallel_2",49.8333339],
#>     PARAMETER["false_easting",150000.013],
#>     PARAMETER["false_northing",5400088.438],
#>     UNIT["metre",1,
#>         AUTHORITY["EPSG","9001"]],
#>     AXIS["Easting",EAST],
#>     AXIS["Northing",NORTH],
#>     AUTHORITY["EPSG","31370"]]

^{Created on 2020-12-03 by the reprex package (v0.3.0)}

Session info

``` r devtools::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.0.3 (2020-10-10) #> os Linux Mint 20 #> system x86_64, linux-gnu #> ui X11 #> language nl_BE:nl #> collate nl_BE.UTF-8 #> ctype nl_BE.UTF-8 #> tz Europe/Brussels #> date 2020-12-03 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2) #> callr 3.5.1 2020-10-13 [1] CRAN (R 4.0.3) #> class 7.3-17 2020-04-26 [4] CRAN (R 4.0.0) #> classInt 0.4-3 2020-04-07 [1] CRAN (R 4.0.2) #> cli 2.2.0 2020-11-20 [1] CRAN (R 4.0.3) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.2) #> DBI 1.1.0 2019-12-15 [1] CRAN (R 4.0.2) #> desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.2) #> devtools 2.3.2 2020-09-18 [1] CRAN (R 4.0.2) #> digest 0.6.27 2020-10-24 [1] CRAN (R 4.0.3) #> dplyr 1.0.2 2020-08-18 [1] CRAN (R 4.0.2) #> e1071 1.7-4 2020-10-14 [1] CRAN (R 4.0.3) #> ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.2) #> evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.2) #> fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.2) #> fs 1.5.0 2020-07-31 [1] CRAN (R 4.0.2) #> generics 0.1.0 2020-10-31 [1] CRAN (R 4.0.3) #> glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2) #> highr 0.8 2019-03-20 [1] CRAN (R 4.0.2) #> htmltools 0.5.0 2020-06-16 [1] CRAN (R 4.0.2) #> KernSmooth 2.23-18 2020-10-29 [4] CRAN (R 4.0.3) #> knitr 1.30 2020-09-22 [1] CRAN (R 4.0.2) #> lattice 0.20-41 2020-04-02 [4] CRAN (R 4.0.0) #> lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.2) #> link2GI 0.4.4 2020-11-30 [1] Github (r-spatial/link2GI@7cbcbdb) #> magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.0.3) #> memoise 1.1.0 2017-04-21 [1] CRAN (R 4.0.2) #> pillar 1.4.7 2020-11-20 [1] CRAN (R 4.0.3) #> pkgbuild 1.1.0 2020-07-13 [1] CRAN (R 4.0.2) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.2) #> pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.2) #> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2) #> processx 3.4.5 2020-11-30 [1] CRAN (R 4.0.3) #> ps 1.4.0 2020-10-07 [1] CRAN (R 4.0.3) #> purrr 0.3.4 2020-04-17 [1] CRAN (R 4.0.2) #> R6 2.5.0 2020-10-28 [1] CRAN (R 4.0.3) #> Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.2) #> remotes 2.2.0 2020-07-21 [1] CRAN (R 4.0.2) #> rgdal 1.5-18 2020-10-13 [1] CRAN (R 4.0.3) #> rgrass7 * 0.2-3 2020-12-03 [1] Github (rsbivand/rgrass7@253f55a) #> rlang 0.4.9 2020-11-26 [1] CRAN (R 4.0.3) #> rmarkdown 2.5 2020-10-21 [1] CRAN (R 4.0.3) #> roxygen2 7.1.1 2020-06-27 [1] CRAN (R 4.0.2) #> rprojroot 2.0.2 2020-11-15 [1] CRAN (R 4.0.3) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.0.2) #> sf 0.9-6 2020-09-13 [1] CRAN (R 4.0.3) #> sp 1.4-4 2020-10-07 [1] CRAN (R 4.0.3) #> stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2) #> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.0.2) #> testthat 3.0.0 2020-10-31 [1] CRAN (R 4.0.3) #> tibble 3.0.4 2020-10-12 [1] CRAN (R 4.0.3) #> tidyselect 1.1.0 2020-05-11 [1] CRAN (R 4.0.2) #> units 0.6-7 2020-06-13 [1] CRAN (R 4.0.2) #> usethis 1.6.3 2020-09-17 [1] CRAN (R 4.0.2) #> vctrs 0.3.5 2020-11-17 [1] CRAN (R 4.0.3) #> withr 2.3.0 2020-09-22 [1] CRAN (R 4.0.2) #> xfun 0.19 2020-10-30 [1] CRAN (R 4.0.3) #> XML * 3.99-0.5 2020-07-23 [1] CRAN (R 4.0.2) #> xml2 1.3.2 2020-04-23 [1] CRAN (R 4.0.2) #> yaml 2.2.1 2020-02-01 [1] CRAN (R 4.0.2) #> #> [1] /home/floris/lib/R/library #> [2] /usr/local/lib/R/site-library #> [3] /usr/lib/R/site-library #> [4] /usr/lib/R/library ```

I did notice some differences however in the WKT string returned by GRASS 7.8.4 and that returned in R (sf 0.9-6), when setting the GRASS projection from a file (with g.proj -c georef=...) / reading the CRS from that file (in sf). I intend to document it in a new issue after a bit more experimentation. UPDATE: probably the rgrass7 repo is not the right place to do that, could be directly GRASS-related.

sf uses GDAL SRS, not PROJ to instantiate WKT. RGDAL uses PROJ rather than GDAL SRS. That may be the reason.

sf uses GDAL SRS, not PROJ to instantiate WKT. RGDAL uses PROJ rather than GDAL SRS. That may be the reason.

I know this is going off-topic, but nevertheless I'l try to summarise my findings on the last point, as it did not turn into a reason to make another issue. The examples refer to the 'Belge 1972 / Belgian Lambert 72' CRS (EPSG 31370).

First of all, in R the WKT returned by sf and by sp are identical:

options(rgdal_show_exportToProj4_warnings = "none")
library(sp)
rgdal::rgdal_extSoftVersion()
#>           GDAL GDAL_with_GEOS           PROJ             sp 
#>        "3.1.3"         "TRUE"        "7.1.1"        "1.4-4"
library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.3, PROJ 7.1.1
wkt(CRS(SRS_string = "EPSG:31370")) == st_crs(31370)$wkt
#> [1] TRUE

^{Created on 2020-12-12 by the reprex package (v0.3.0)}

This WKT string (WKT2) is identical to the one returned (in the shell) by gdalsrsinfo EPSG:31370 -o wkt2_2018 and by projinfo EPSG:31370 -o WKT2:2019
Same equality found between gdalsrsinfo EPSG:31370 -o wkt1 vs. projinfo EPSG:31370 -o WKT1:GDAL (WKT1) and between gdalsrsinfo EPSG:31370 -o wkt_esri vs. projinfo EPSG:31370 -o WKT1:ESRI (ESRI WKT, although printed as a one-line string by projinfo).
GRASS 7.8.4 g.proj -w epsg=31370 returns WKT1, almost exactly the same as the WKT1 from gdalsrsinfo / projinfo, the only difference is that the g.proj output has an extra TOWGS84 line. I assume WKT2 is expected in future GRASS versions. The WKT1 in GRASS vs WKT2 in R contrast does not pose problems when exchanging data objects with sp/sf packages in R, which is great.

But I did see a clash between objects with the 'Belge 1972 / Belgian Lambert 72' CRS. However it solely had to do with a WKT difference between a location created in GRASS as EPSG:31370 and official spatial data sources in Flanders: the latter often contain ESRI's take on 'Belge 1972 / Belgian Lambert 72' (referred on a reference list by ESRI as 'Well-known ID 31370'). Although it may differ negligeably from EPSG:31370 with respect to coordinate accuracy, it does differ syntactically and numerically: switch of the standard parallels and different numerical precision of false Easting/Northing. I don't know how common such kind of ESRI vs EPSG difference is. My conclusion is that all I can do (locally) is overriding ESRI's CRS with EPSG:31370.

If you think this needs rediscussion elsewhere, I'll be glad to post it. Just I didn't feel it was worth it.

With GRASS 7.8.5dev built against PROJ 7.2.0 and GDAL 3.2.0 (on PROJ 7.2.0), the output of g.proj -w epsg=31370 matches PROJ and GDAL WKT1, lines2:24. Differences may come from PROJ versions. No TOWGS84 should be shown in modern PROJ, as transformation pipelines should be used. g.proj should arguably emit WKT2:2019. From https://github.com/OSGeo/grass/blob/master/general/g.proj/output.c, I think that some code is present, but I cannot see how to request a format in g.proj.

The file you refer to was indeed changed in the course of PR OSGeo/grass#976; these changes are present only in the master branch. In this PR, it is stated (https://github.com/OSGeo/grass/pull/976#issuecomment-694491675):

WKT is the latest WKT2 version if compiled against GDAL3+, otherwise WKT1.

The last change of this file (on master) was on 24 Sep.

But it is not part of the branch with the 7.8 releases, where the latest change of the file appears to be from 8 Mar 2019. So that's probably the reason why it returns WKT1.

PR 976 was added to the 8.0.0 milestone, which I assume is what the master branch is going to become.

rsbivand / rgrass

gmeta() returns a proj4string #15