How to define a local equidistant CRS
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 at nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en