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

[Roger Bivand]

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

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. The package is not on 
CRAN, so was not caught by our reverse dependency checks, so the 
responsibility for making sure things work rests on the maintainer, not 
us.

Please also convert the example to points, (you reach 
.spTransform_Polygon() here, so these are not points), and try proj from 
the command line. I think that your +a and +b also have units problems (km 
not m). The +R_A argument is not given as in current:

https://proj.org/operations/projections/aeqd.html

> In addition: Warning messages:
> 1: In showSRID(uprojargs, format = "PROJ", multiline = "NO",
> prefer_proj = prefer_proj) :
>  Discarded ellps unknown in CRS definition: +proj=aeqd +a=6378.137
> +b=6356.752 +R_A +lat_0=30.4 +lon_0=35.1 +x_0=0 +y_0=0 +type=crs
> 2: In showSRID(uprojargs, format = "PROJ", multiline = "NO",
> prefer_proj = prefer_proj) :
>  Discarded datum unknown in CRS definition
>
> If I change the CRS definition as follows:
> projstring <- paste("+proj=aeqd +R_A +lat_0=", YMiddle,
>                    " +lon_0=", XMiddle,
>                      " +x_0=0 +y_0=0 +ellps=WGS84",sep="")
>
> Then using proj 6.3.2 we get a warning:
> "In showSRID(uprojargs, format = "PROJ", multiline = "NO", ... :
> Discarded datum Unknown based on WGS84 ellipsoid in CRS definition"
>
> whereas on a newer setup, with "Loaded PROJ runtime: Rel. 7.2.1" the
> above definition is accepted silently.
>
> Our question is what is the reliable way to define a local Cartesian
> CRS which would work anywhere in the world, and preferably be
> equidistant.

With a reproducible example that works with proj command line tools, one 
could ask for views on the PROJ list itself.

Hope this helps,

Roger

> Should we just find the local UTM zone and use that?  UTM
> is not equidistant, but it is ubiquitous, and if the cellular network
> is not too broad, I assume that distances will be fairly accurate.
>
> Best regards,
> Micha
>
> [1] https://github.com/overeem11/RAINLINK
>
>

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