Tutorial(s) and communicating on using PROJ>=6/GDAL3 in R #43
Comments
|
Thanks, useful. I'll be speaking at whyR 2020 on Sunday 27 September at 11.00 - some parts of the talk and the (to be posted soon) code extend the celebRation 2020 talk you reference. I'll post a link once I have one. I haven't read through your tutorial carefully yet, but a question that is becoming more visible is that EPSG:4326 is by definition Latitude-Longitude, while OGC:CRS84 is by definition Longitude-Latitude. From the GDAL barn-raising, PROJ/GDAL are trying to be authority-compliant, but this isn't resolved yet on our side. In the talk, I pivot to using OGC:CRS84 instead of EPSG:4326; you may think about what works for your target readership, but probably EPSG:4326 is a bad idea for GIS/visualization from the point of view of axis order. |
|
That's very interesting to learn about @rsbivand , it was beyond my knowledge (I thought we could always reside to EPSG). I will definitely have look at that - looking forward to your talk. |
|
I agree that writing |
|
Link to my Sunday whyR talk (11.00 CEST) probably https://youtu.be/mEAyQ8bv1zU. Followed by interesting panel. Presentation files on https://github.com/rsbivand/whyR20_files; output from current rsbivand/sp and rev 1066 of rgdal, sf 0.9-6 CRAN. |
|
Thanks Roger for the interesting extra slides (42-73), which you couldn't discuss in the talk, about (among others) proj4string stripping, axis order, WKT representation by GDAL versus PROJ (and, some new |
|
It is helpful, except, I think, for the statement that one should "specify the CRS by using the EPSG code, but do so without using a proj4string". That is not generally usable advice. While you can use these codes for the few established CRSs, you cannot use it to define an appropriate system for the data you are working with (and there is an infinite number of these). While we can still use PROJ notation, at least with the WGS84 datum, that seems to be the only generally usable approach available; unless you are able to create your own WKT. I am not aware of any effort to write software to support creation of WKT from basic components (like a PROJ string), but that could be very useful. Either way, it would be good if your tutorial explained how to do that, if there is a way. |
indeed, in fact I've been considering writing out a short explanation of centring a local projection on a lon,lat location and using metres to buffer an extent around it. It's handy to show the cosine hack one would use to get an approximately correct aspect ratio with angular coordinates, but a throw way proj-string projection and a metres-either-side range is much simpler. (sometimes you care more about arbitrary centre and extent than you do about the projection family, and in addition you may not care at all about epoch-related inaccuracy). I've been using this as a target spec with the GDAL warper and it's clear it's the best way to generate a local projection for general use. It could be expanded into a "we use this epoch in such-and-such context", so here's the WKT2 version of a throw away proj string used in the more casual ways .... |
|
I may be wrong but think the issue is the difference between using proj4 strings as coordinate reference systems, or as coordinate transformations. If we use them as CRS, they get |
|
Yes, except code is also representation, and round-tripping has no hope to exist. The problems are real ones, and the benefits are clear only in local and primarily terrestrial contexts, exactly the ones where a local authority grid likely already exists and exploring the properties of projections has less value. The WKT2 rigour in the closed contexts of sf and terra with their built-in GDAL stacks is unproblematic afaic, they just aren't the entire story. |
If you set up your own CRS parameters, changes are that this is what you want.
But what other option would I have? That way I could get at least get a WKT2 for the WGS84 variant; and then I could try to use that as template by replacing WGS84 with the datum of choice. But I have never needed that; for the reasons that Mike alludes to.
To extend this discussion a bit: when using the WGS84 datum, the So is it fair to say that:
? |
|
That summary accords with my understanding, which has been slow - you can also use the modern equivalent of I think the risks are at sub-metre level for ensemble datums, if you don't have the time stamp of your coordinate measures and ability to choose the right epoch. Other risks are the standard ones we always had (100s or so metres for failure to specify the actual datum). You also can't use |
|
Is it possible that an R function might be needed for custom CRS to either just write a WKT2-2019 or PROJJSON string from user Proj4-like input? If PROJ was available, the representation could be validated, and/or I recall that PROJJSON was suggested as somewhat easier to handle than WKT2, but for user input, Proj4 strings are clearly easiest. The downside is of course that EPSG coded CRS are much easier to fit into SQL search for transformation pipelines. I've seen 1-2m given as the EPSG:4326 ensemble inaccuracy, with further comments on https://lists.osgeo.org/pipermail/gdal-dev/2021-May/054112.html , https://github.com/rouault/gdal/blob/rfc_81/gdal/doc/source/development/rfc/rfc81_coordinate_epoch.rst . I don't think things have settled at all. If we were already an OSGeo community, maybe we'd be able to put forward points of view, I'm unsure about that, though. |
|
Thanks @rhijmans for the feedback! I'd be happy to improve the tutorial where possible. Related to @mdsumner 's remark, its scope is indeed limited to officially registered (say EPSG) CRSs - probably we should state that more clearly. But I feel it is the scope most users are confronted with. So importantly, the tutorial was not intended to explain on defining CRSs from scratch. I think that would best fit in a separate tutorial - cf. post above. The tutorial's aim is to give generally applicable, basic advice, easy for (beginning) R users (in the light of the recent proj4string obsoletion). In that regard (keeping things simple and robust), I doubt whether it's a good idea to showcase that one can also use IMO the easiest / failsafe way in specifying an 'official' CRS (for spatial data) is referring the CRS from an official database (like
Thanks! In order to make the tutorial more complete, it may indeed be a good idea to add a short note on the WGS84 datum case in parentheses or as a footnote (referring to
Sure one could make mistakes and pass the wrong key. But I think it's a better managed approach - to comply with EPSG CRSs that is. I don't understand 'harder to read' in your explanation.
There is much merit indeed in explaining more on what can be done (and not) with projstrings. It is not an area I'm familiar with though. IMO this deserves its own tutorial (you're also welcome to add it to the inbo/tutorials repo, even though it's linked to my organization, but it would then sit next to the EPSG-CRS tutorial). -- As a sidenote, and referring to aforementioned WGS84 inaccuracies, IMHO a CRS with a specific WGS84 datum realization (latest realization is G1762: EPSG:1156), such as the EPSG:9057 CRS, seems preferable (when possible) to use over EPSG:4326 (or OGC:CRS84) CRS which uses the WGS84 ensemble. The latter (datum:EPSG::6326) is an ensemble datum (PS: I'm not often online these weeks, sorry for the late reply) |
|
My main point was that saying that the general principle is to
seems extreme and exaggerated. Instead you could recommend to use EPSG codes when available; and note that you can use PROJ.4 notation as long as you do not stray from WGS84. There is, by the way, nothing "official" about these codes. I am not suggesting to explain how to roll your own WKT2 CRS with another datum (as we do not seem to know how to easily do this!)
There might be many that have other needs. I have worked a lot with species distribution data. These do not tend to follow national boundaries and so the appropriate CRS is not in the databases. There is a national grid for Belgium, one for the Netherlands and (presumably) one for Luxembourg. What do you use when you map the Benelux? So I think the situation in which you want to define your own parameters is very common. Of course you could pick one that is "close enough" from a list, but is that good advice?
To me this is hard to read: ((CRSs create a lot of misunderstanding. Using codes does not educate. It would be have been much better to expand PROJ notation, it seems to me. I do not know the limits of that, but the point is moot, we are going to have to live with it.))
How would that matter in data analysis (as opposed to the e.g. the perfect drone delivery service, be it of goods or bads); as long as the original datum is known and the data correctly transformed to the same (appropriate) datum? |
|
Yes, this is so important. The idea that you must find a code is incredibly toxic, and you can see this idea propagated online in so many communities where experienced users should know better (it's similar to the "find your UTM zone" advice, which is hopelessly naive and off-base). |
|
If it the case that the position data are only going to be used by the person choosing the CRS using a custom Proj4 string, without any datum transformation, then by all means use Proj4 strings. If the data are going to be written to file/db, then the CRS representation matters much more, as is the case for datum transformation. Another question - for vector data, an s2 approach with only GEOGCRS on a sphere is probably easiest - then everybody is in the same situation. But how do raster and s2 work together? This is, I feel, somewhere in the gdistance area, isn't it? s2 to WebMercator isn't accurate, but at scales of several hundred km, does it matter? Potentially yes? If so (data integration from multiple sources), we are still stuck. Admittedly, GRASS locations could be defined somewhat like Proj4, but in that case, one transformed all non-compliant input to the location CRS. Do we still need to do that, or is s2 spherical without worrying about epochs or datum ensembles close enough? I'd be very interested in seeing @paleolimbot 's views. For species distribution data, is a planar representation essential, or is a spherical representation adequate, and side-steps the CRS representation question? |
|
Honoured that anybody is interested in my views! s2 shines when the data are already provided as longitude/latitude, which is and I imagine will be common for a long time. Correcting for epoch or other CRS issues will probably always be at the beginning or end of a script...s2 will always sit in the middle. The other place where s2 shines, regardless of the initial form of the data, is when you need to ask a global question (e.g., "within 500 km of a coastline but not on land"). Those questions rarely have an accuracy requirement that is impacted by the spherical assumption. Maybe a good way to summarize my feelings about this are that there are no "innaccurate" reference systems, only inaccurate transformations, inaccurate distances between points, and inaccurate representations of edges (i.e., how you draw a line between two coordinates). s2 can help by minimizing the number of times one has to invoke a transformation and has tools that can make edge representation explicit by adding points (at the expense of all distances in s2 being innaccurate/spherical only). |
|
@paleolimbot thanks! My feeling is that maybe we can lessen the need for custom projections by choosing to use spherical representations to a greater extent from now on, rather than having to go planar to handle distances/areas/intersections. However, perhaps species distribution models, if they presume planar, might need re-visiting? Are perhaps also global spherical grid systems, like the now-archived dggridR package https://github.com/r-barnes/dggridR/ or similar? I feel that, unless analysts need high precision, spherical is not worse than planar some distance from the central lon-lat point. |
|
I agree that using lon/lat can be very good for some types of data analysis, but that is a different discussion. There are many instances in which someone would reasonably want to define a custom CRS. As-is, in the GDAL/PROJ world we only have the proj4 string for that, even if its usability has degraded. In the tutorial, Floris says that the proj4 notation should not be used, in part because he believes it will likely be removed. That may be true or not, but if we could only use EPSG (etc) codes we would not have competent general purpose spatial data analysis software. We would need additional software that translates a (perhaps expanded) proj4 string (or similar) to wkt2. (and, who knows, it could then be included in the next version of PROJ? ) |
|
I think all these cases are already covered in both sf and sp workflows. The relevant problem is that GDAL degrades input file CRS when exporting to Proj4, so the internal R-side CRS representation in sp or sf workflows cannot be Proj4 when reading from file or checking the representation for validity against PROJ. PROJ and GDAL both import Proj4 strings, but deprecate with default insertion of WGS84 as datum in the WKT2-2019 representation if So R-side WKT2-2019 can be created from Proj4 string user input, usually without an authority for the PROJCRS. The WKT2-2019 should also be OK to write to file as part of an object, and should be retrieved as written. It would be useful to check transformation pipelines between different |
|
My take, merely for transformations available to general purpose software as Robert mentioned, is to have a basic interface to what PROJ does. @paleolimbot provides this in libproj but has cryptic and painful problems keeping that on CRAN, they have been difficult and unhelpful with failed communications from their side and failed to indicate an impending archive step because of problems in the package. Communications have been ambiguous, and it depends who on the CRAN team is replying. That's typical in my experience, the more in the club you are the better, with no transparency about what they support, help with, etc. They try to be robotic, objective, but they are not. If you can use PROJ as-is, you get everything you need, input proj4, wkt, whatever format the library provides, and the user/developer can make choices. If we had that shared resource then it'd be easy to write an adjunct to the blog post discussed. |
|
pkg libproj made the proj sources part of the package, which sounds attractive because you loose the need to install proj besides R when doing source installs, but which also puts all the PROJ sources under the CRAN src checking, leading to problems similar to e.g. those we now face with s2 (g++11 issues in the upstream library, us needing to push google to move for our case, everyone waiting). I don't see why the way libproj interfaced PROJ can't be done without including all of PROJ in src but simply using the way sf/rgdal/terra etc use PROJ. If that is what you need then why don't you do that? I also don't see how complaining about CRAN contributes to the discussion here. |
|
This also relates for PROJ to the size of its distributed data; having a projdata package has been suggested, but is version-fragile as Windows and MacOS CRAN binaries may use different versions. This also impacts OSGeo4W, as users may see one PROJ version in R (even if looking at the OSGeo4W-bundled data) and a different version of the data. This matters now because EPSG 10 changed the definitions of the database table and field structures in some places (so R/C code built against later PROJ through PROJ itself bites the dust when looking at an earlier DB structure). |
|
Thanks @rhijmans for shedding light on your usecase (or data analysis in general), and for advocating being able to define one's own CRS (I like the Benelux example 🙂). Currently I'm unsure what's the best advice to easily define a new CRS (see below) and this may need discussion at PROJ directly. It seems it is not currently a purpose of PROJ to make this easy. Enlisting here some things noted, after reading through pages on https://proj.org (but nothing more):
All in all, it seems that the particular case of creating a new CRS is not specifically envisaged, but I may be wrong. The general idea seems to be that an actual CRSs should be defined in full (WKT or equivalent). So this ultimately seems a debate between full versus partial CRS specifications.
Anyway, I intend to add a bit of nuance to the tutorial, regarding projstrings. Probably it's safer to say 'goodbye PROJ.4' than 'goodbye proj4strings' 😉. |
|
Lots of good stuff here! All my bits are off topic...I'll respond here but feel free to open threads elsewhere if there's any interest in further discussion.
|
|
After reading further in PROJ documentation I conclude that some thoughts in my previous post, about PROJ assuming an implicit datum, and hence confusion in the documentation, were incorrect - sorry for that. I'll try to clarify that part here (corresponding parts are crossed out above). As I understand now, (modern) proj-strings only define the coordinate operation itself, i.e. without defining the input and output CRS (and hence, their datum). So if a datum shift is needed (because input & output CRS have different datum), then the corresponding shift must be defined in the pipeline, but the actual CRS datums don't. Another way to say this would be that the geoid's position relative to the ellipsoidal CS are (generally) neglected in proj-strings. For operations with geographical coordinates only the ellipsoid model and relative positions between ellipsoidal CSs (e.g. by Helmert transformation) are defined by proj-strings: one doesn't need the geoid position (relative to the ellipsoidal CS) to define these coordinate operations. Of course you still do need the datum of your input and output CRS in order to match coordinates with real locations. If I remember correctly, PROJ.4 did assume a WGS84 datum default and allowed an alternative datum definition (e.g. directly with Meanwhile, I'm preparing some questions, to address at PROJ, about the role of proj-strings in declaring CRSs (including custom ones). |
|
Below follow some results from attempts to declare a CRS from a proj-string, with recent PROJ (using -- While the declaration of a CRS using a proj-string seems not documented on proj.org for the main PROJ programs (
(*) 'little support for datum': that is, except when you're happy to use a BOUNDCRS, which is generated when adding a
IMO the last steps (support for BOUNDCRS with For the R side I understand that "The relevant problem is that GDAL degrades input file CRS when exporting to Proj4" (Roger's comment) so this would mean that adding |
|
Recent sp/rgdal provide where the C code is: BOUNDCRS are found most often when a file contains a |
|
Thanks for pointing at the argument. Since for below example the argument does not have an effect, it seems as if the Reprexlibrary(sp)
rgdal::rgdal_extSoftVersion()
#> GDAL GDAL_with_GEOS PROJ sp
#> "3.2.1" "TRUE" "7.2.1" "1.4-5"
# old proj4string (with towgs84) for EPSG:31370:
##################################################
crs1 <- CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs +type=crs",
get_source_if_boundcrs = TRUE)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on International 1909 (Hayford)
#> ellipsoid in Proj4 definition
# same, but get_source_if_boundcrs = FALSE
##################################################
crs2 <- CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs +type=crs",
get_source_if_boundcrs = FALSE)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on International 1909 (Hayford)
#> ellipsoid in Proj4 definition
# current proj-string (without towgs84) for EPSG:31370:
##################################################
crs3 <- CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m +no_defs +type=crs",
get_source_if_boundcrs = TRUE)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on International 1909 (Hayford)
#> ellipsoid in Proj4 definition
# same, but get_source_if_boundcrs = FALSE
##################################################
crs4 <- CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +units=m +no_defs +type=crs",
get_source_if_boundcrs = FALSE)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
#> = prefer_proj): Discarded datum Unknown based on International 1909 (Hayford)
#> ellipsoid in Proj4 definition
all.equal(crs1, crs2)
#> [1] TRUE
all.equal(crs1, crs3)
#> [1] TRUE
all.equal(crs1, crs4)
#> [1] TRUE
cat(wkt(crs1), sep = "\n")
#> PROJCRS["unknown",
#> BASEGEOGCRS["unknown",
#> DATUM["Unknown based on International 1909 (Hayford) ellipsoid",
#> ELLIPSOID["International 1909 (Hayford)",6378388,297,
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]]],
#> CONVERSION["unknown",
#> METHOD["Lambert Conic Conformal (2SP)",
#> ID["EPSG",9802]],
#> PARAMETER["Latitude of false origin",90,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8821]],
#> PARAMETER["Longitude of false origin",4.36748666666667,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8822]],
#> PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8823]],
#> PARAMETER["Latitude of 2nd standard parallel",49.8333339,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8824]],
#> PARAMETER["Easting at false origin",150000.013,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8826]],
#> PARAMETER["Northing at false origin",5400088.438,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8827]]],
#> CS[Cartesian,2],
#> AXIS["(E)",east,
#> ORDER[1],
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]],
#> AXIS["(N)",north,
#> ORDER[2],
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]]]Created on 2021-06-18 by the reprex package (v2.0.0) Session infosessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.1.0 (2021-05-18)
#> 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 2021-06-18
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date lib source
#> cli 2.5.0 2021-04-26 [1] CRAN (R 4.1.0)
#> digest 0.6.27 2020-10-24 [1] CRAN (R 4.1.0)
#> evaluate 0.14 2019-05-28 [1] CRAN (R 4.1.0)
#> fs 1.5.0 2020-07-31 [1] CRAN (R 4.1.0)
#> glue 1.4.2 2020-08-27 [1] CRAN (R 4.1.0)
#> highr 0.9 2021-04-16 [1] CRAN (R 4.1.0)
#> htmltools 0.5.1.1 2021-01-22 [1] CRAN (R 4.1.0)
#> knitr 1.33 2021-04-24 [1] CRAN (R 4.1.0)
#> lattice 0.20-44 2021-05-02 [4] CRAN (R 4.1.0)
#> magrittr 2.0.1 2020-11-17 [1] CRAN (R 4.1.0)
#> reprex 2.0.0 2021-04-02 [1] CRAN (R 4.1.0)
#> rgdal 1.5-23 2021-02-03 [1] CRAN (R 4.1.0)
#> rlang 0.4.11 2021-04-30 [1] CRAN (R 4.1.0)
#> rmarkdown 2.8 2021-05-07 [1] CRAN (R 4.1.0)
#> rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.1.0)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.1.0)
#> sp * 1.4-5 2021-01-10 [1] CRAN (R 4.1.0)
#> stringi 1.6.2 2021-05-17 [1] CRAN (R 4.1.0)
#> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.1.0)
#> withr 2.4.2 2021-04-18 [1] CRAN (R 4.1.0)
#> xfun 0.23 2021-05-15 [1] CRAN (R 4.1.0)
#> yaml 2.2.1 2020-02-01 [1] CRAN (R 4.1.0)
#>
#> [1] /home/floris/lib/R/library
#> [2] /usr/local/lib/R/site-library
#> [3] /usr/lib/R/site-library
#> [4] /usr/lib/R/library |
|
Interesting. I cannot check without re-installing sp and rgdal from CRAN, as the development versions (''sp** my github fork, rgdal R-forge) do different things, including caching CRS objects (some packages make repeated calls to |
|
I don't see any difference after re-installing sp and rgdal from CRAN. The The chronology by reconstruction was possibly that I switched all instantiation via |
|
In development sp, identical results are caused not by a logic error, but by (new) CRS look-up caching. Repeated look-ups retrieve output CRS from a local cache if the same input is received. |
|
Thanks a lot Roger for demonstrating the effect of In a context of changing use of proj-strings, IMHO I also think it is natural to consider PROJ as reference. Currently it still supports Some things to note:
|
|
Sharing here some outcomes from the PROJ mailing list, earlier this week, on the subject of using PROJ strings to specify a CRS (answers from Even Rouault and others). PROJ strings for CRSs will get continued (though reduced) support, for backward compatibility reasons. After clarifying some aspects, following advice is considered appropriate (see also https://lists.osgeo.org/pipermail/proj/2021-June/010301.html):
I intend to use that as basis for a few updates in our tutorial on CRSs in R. Further notes:
|
Adding some insights on the current and expected state of PROJ string support, since several people stick to PROJ strings for convenience despite this approach is discouraged. See r-spatial/discuss#43 (comment) and https://lists.osgeo.org/pipermail/proj/2021-June/010301.html .
|
Just updated the discussed CRS tutorial (link). I essentially added a bit more info on current state of PROJ string support for CRSs (partly in footnotes), and referred to PROJJSON, meanwhile trying to stay in line with PROJ documentation and trying not to complicate things too much. Basic advice for specifying EPSG CRSs of course remained the same. |
|
Maybe replace Also consider introducing something like: to help in finding EPSG codes. Another optional extra is: The latter needs updating to link directly to https://www.asprs.org/wp-content/uploads/2018/03/10-16-GD-Belgium.pdf for the 2016 link from https://www.asprs.org/asprs-publications/grids-and-datums. These short articles give useful context on how countries have adapted and modernised their standards. |
OK, thank you Roger.
Yes, that's a good fit indeed! The ASPRS Grids and Datums articles are an interesting source, thanks for reminding. IIRC for Belgium it focuses less (resp. not) on the two currently used CRSs: 1) "Belge 1972 / Belgian Lambert 72" (EPSG:31370, based on the "Reseau National Belge 1972" datum) - still used most often - and 2) "ETRS89 / Belgian Lambert 2008" (EPSG:3812, based on the ETRS89 datum) which has accuracy advantages in using GPS data and for conversion to other ETRS89-based CRSs within Europe (EPSG or custom CRSs). In time, it can be expected that the second will have more priority; it is also required by the EU INSPIRE directive. (Comments by the Flemish Geographic Information Agency are at https://www.geopunt.be/voor-experts/lambert-2008 - only available in Dutch). I hope to delve deeper into this later on (consulting official Belgian sources, including related topics such as the Flemish Positioning Service FLEPOS and the unfortunate ESRI:102499 mismatch), and perhaps start a separate page about it, for Belgian R users and other interested people 😉. |
Inspired on: - a (long-standing) suggestion of Roger Bivand in r-spatial/discuss#43 (comment) - since rgdal will retire, the hint at the end of an email by Michael Sumner: https://stat.ethz.ch/pipermail/r-sig-geo/2021-September/028778.html The rgdal approach has been included as a markdown comment.
This applies the request by Roger Bivand in r-spatial/discuss#43 (comment). Apologies for this late update.
Opening a new discussion, as I've been trying to write a tutorial on current good practice in specifying a CRS in R (with focus on the commonly used
sf,spandraster). It would be great if some people could have a critical look at this tutorial. Comments / pull requests are most welcome, and that includes language improvement. Getting feedback on it was the primary reason for me to open this issue.But, feel free to widen up discussion, about communication (needs?) to R users on practical PROJ>=6/GDAL3 aspects they must know about. (And, who knows, I may have missed existing communications so I'd be happy to hear about it!) While currently most developers and more advanced R users will be aware of these changes, and have found where to look on the internet, I think that many R users are still unaware (even though they should be alarmed by the proj4strings warnings). At least that's my experience - after all these are recent changes.
The context of the above linked tutorial is the following. In our organization (Research Institute for Nature and Forest - www.inbo.be/en) we are gradually extending a tutorials website (https://inbo.github.io/tutorials/), mostly oriented towards reproducible science (spatial topics are found by clicking on the 'gis' tag.). Its primary aim is to support and stimulate scientists within the organization, but of course it is available to everyone and open to contributions. And next to the website, some very nice colleagues organize coding clubs in order to accelerate hands-on experience with R. Overall, we have about 150 R users (as derived from our internal useR mailing list), with a large variety in their R experience and in their effective R usage. Not all scientists are enough self-reliant to find out about technical aspects they need to know, hence they desire being supported in some coordinated way. I imagine this is within normal expectations and it is the case elsewhere.
In order to fill that gap we believe organization-level tutorials are of help. In some cases we just make a tutorial to provide a link to another tutorial elsewhere, in order to prevent unnecessary duplication and keep maintenance low.
The text was updated successfully, but these errors were encountered: