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

[Roger Bivand]

On Fri, 22 Jan 2021, 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.
>>>>>    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):

Could a longer-term solution be to  use the s2 package to do all of your 
operations on the sphere? Then you avoid having to find a local 
projection.

Roger



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

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