[postgis-users] Implementation of postgis interoperability with R

[Etienne B. Racine]

Dear list,

I'm currently updating the simple feature <https://github.com/r-spatial/sf>
package for R to use the latest proj6 to interact with postgis. One of the
things we perform behind the scene is synching the projections between R
and postgis.

The current version of the package provides interoperability in two ways:
read and write data to postgis, and allow the user to write R code that is
executed, as SQL, on the remote database, just like a normal data.frame (R
tables).

Our current implementation relies on epsg code and proj4string matching to
synch database srid with the local crs, but with PROJ6, we need to add proj
wkt matching.

Two questions:

   - This is how I plan to adapt the projection synching, any suggestion
   for improvement or pointer to other implementations would be much
   appreciated:
   1. Match local epsg code and srid, create a crs in R using the `srtext`
   and check that the two projections are equivalent using `IsSame` from
   OGRSpatialReference.
   2. No matching code, so check `spatial_ref_sys` table to see if a proj
   WKT would match any `srtext` (slow and not very robust, I have a follow up
   question)
   3. No matching code and no WKT match, so try to match proj4string.
   4. Clearly no matching projection, so insert new srid in
   `spatial_ref_sys` and set `srid = max(srid) + 1`, `auth_name = 'sf'` (the
   simple feature package in R is called sf).
   5. If it fails (e.g. lack of permission), then error and ask the user to
   change the projection, or use srid = 0.


   - Is string matching `=` the best way to match a `srtext` and
   `proj4string` in  `spatial_ref_sys`?
   ``` sql
   select * from spatial_ref_sys where srtext = {wkt};
   ```
   But `=` obviously rejects matches if parameters are ordered differently
   (and I suspect we could have issues with encoding). One solution is to read
   `spatial_ref_sys` to R, and match there, but it is very very slow since it
   requires to parse all the wkt to create a projection, and then perform a
   comparison with OGRSpatialReference->IsSame on every projection.

Thanks for your help!
Etienne
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20200312/f952626a/attachment.html>