How to define a local equidistant CRS

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.
   Among the functions in RAINLINK, some require to get the attenuation
   for each microwave link (pair of towers), from those towers that are
   nearby, within, say 15 km. This means that the tower locations need
   to
   be transformed to an equidistant Cartesian CRS.

   In the current version of RAINLINK, this was done as follows:
   An sp object was prepared from the tower locations, with CRS
   "EPSG:4326". The average latitude and longitude values were obtained
   as "YMiddle" and "XMiddle". Then a new, one-off azimuthal equidistant
   CRS was defined as:
   # Set projection string
   projstring <- paste("+proj=aeqd +a=6378.137 +b=6356.752 +R_A
   +lat_0=",
   ?? YMiddle, " +lon_0=",
   ?? XMiddle," +x_0=0 +y_0=0",sep="")

   This worked until the recent advances in proj and gdal. Now with proj

   =6.x and gdal >=3.x the above definition causes an error:

   Error in .spTransform_Polygon(input[[i]], to_args = to_args,
   from_args
   = from_args,? :
   ?coordinate operation creation from WKT failed: generic error of
   unknown
   ?origin

  ##=====================================

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:

 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

SpatialPoints:
          x        y
[1,] 6.98795 53.09842
[2,] 6.94800 53.08685
Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
+no_defs

 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))

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:

 options("rgdal_show_exportToProj4_warnings"="none")
 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

SpatialPoints:
          x        y
[1,] 6.98795 53.09842
[2,] 6.94800 53.08685
Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
+no_defs

 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))

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

at least until sp 1.5 and rgdal 1.6, when the default will change to suppress 
these warnings.

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))


If this isn't a reprex, please extend briefly.

Roger