|
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.
If this isn't a reprex, please extend briefly.
Roger
Thanks again,
Regards, Micha
--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
|
