[R-sig-Geo] How to define a local equidistant CRS

[Roger Bivand]

On Sat, 23 Jan 2021, Micha Silver wrote:

> 
> Roger:
> 
> 
> Thanks for your double effort to help.
> 
> Sorry about the line breaks.
> 
> (please see inline)
> 
> 
> On 1/22/21 6:19 PM, Roger Bivand wrote:
>       On Fri, 22 Jan 2021, Roger Bivand wrote:
>
>             On Fri, 22 Jan 2021, Micha Silver wrote:
>
>                    Hello Roger:
> 
>
>                    On 1/20/21 12:39 PM, Roger Bivand
>                   wrote:
>                           On Wed, 20 Jan 2021, Micha
>                         Silver wrote:
>
>                                 The RAINLINK
>                               package [1]
>                               derives rain
>                               rate from the
>                               attenuation of
>                                 signal
>                               strength between
>                               microwave
>                               towers. Data
>                               typically comes
>                               from
>                                 cellular
>                               network
>                               providers, and
>                               the location of
>                               the towers is
>                                 typically
>                               given in
>                               longitude/latitude. 
>
>                           This is not a reproducible
>                         example. Please provide a
>                         reproducible
>                           example,
>                           specifying all of the
>                         software versions,
>                         especially PROJ and GDAL.
>                         The
>                           changes in PROJ are now
>                         entering the 4th year, so
>                         I'm a bit surprised
>                           that
>                           steps were not taken
>                         before to check for problems
>                         to your workflow as
>                           Proj4 strings shift from
>                         deprecated to defunct.
> 
> 
> 
> 
> This brings us back to the original question:
> 
> 
> What would be the recommended, reliable way to create an equidistant CRS on
> the fly, for any location determined by a cluster of points in Long/Lat??
> and avoiding using proj.4 strings??
> 
> Should we:
> 
>     - construct a WKT description to create the CRS?
> 
>     - just revert to using UTM (not equidistant, but available everywhere,
> and easy to determine the zone, EPSG code for any location)
> 
>     - Something else?
> 
> 
>
>                    
>
>       If this is a minimal reprex (making "towers" from the first row,
>       and putting +lat_0= and +lon_0= at reasonable values):
>
>       library(sp)
>       df <- data.frame(x=c(6.98795, 6.948), y=c(53.09842, 53.08685))
>       coordinates(df) <- c("x", "y")
>       slot(df, "proj4string") <- CRS("EPSG:4326")
>       df
>       projstring_aeqd <- paste("+proj=aeqd +lat_0=53.09 +lon_0=6.97",
>        "+x_0=0 +y_0=0 +ellps=WGS84")
>       spTransform(df, CRS(projstring_aeqd))
>
>       then I don't see any problem:
>
>       SpatialPoints:
>                    x         y
>       [1,]  1202.371  937.1959
>       [2,] -1474.053 -350.3307
>       Coordinate Reference System (CRS) arguments: +proj=aeqd
>       +lat_0=53.09
>       +lon_0=6.97 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
>       Warning message:
>       In showSRID(uprojargs, format = "PROJ", multiline = "NO",
>       prefer_proj = prefer_proj) :
>         Discarded datum Unknown based on WGS84 ellipsoid in Proj4
>       definition
>
>       The warning just says that you need to know about and decide how
>       to handle changes in CRS representations. Maybe you do not
>       realise that sp::CRS() and sp::spTransform() both call rgdal if
>       available, so you can use:
>
>       options("rgdal_show_exportToProj4_warnings"="none")
>
>       before loading sp to mute warnings:
> 
> 
> 
> OK, clear
> 
>
>       Alternatively, use the equivalent:
>
>       library(sp)
>       df <- data.frame(x=c(6.98795, 6.948), y=c(53.09842, 53.08685))
>       coordinates(df) <- c("x", "y")
>       slot(df, "proj4string") <- CRS("EPSG:4326")
>       df
>       projstring_aeqd <- paste("+proj=aeqd +lat_0=53.09 +lon_0=6.972",
>        "+x_0=0 +y_0=0 +datum=WGS84")
>       spTransform(df, CRS(projstring_aeqd))
> 
> 
> So specifying the "+datum=" in the proj.4 string also avoids the warning.
>

Yes, but the exceptions for WGS84, NAD83 and NAD27 are just that, 
exceptions from the general rule of +datum= being dropped. The next road 
bump will occcur soon(ish) because WGS84 is an ensemble of datum 
specifications, not a unique specification, and PROJ 7.2.0 is starting to 
address this: https://proj.org/news.html (Update to EPSG 10.003 and make 
code base robust to dealing with WKT CRS with DatumEnsemble (#2370)). The 
link to https://github.com/OSGeo/PROJ/pull/2370 gives details.

Unless you need sub-metre accuracy, this may not matter, but if different 
parts of the world use different WGS84 definitions (at different times), 
the s2 package may be more stable though a little less accurate.

This is more detail than you expected, I think, but at least it is a way 
of documenting things.

Roger

> 
>
>       If this isn't a reprex, please extend briefly.
>
>       Roger
> 
> 
> Thanks again,
> 
> Regards, Micha
> 
> 
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en