[R-sig-Geo] How to define a local equidistant CRS
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Fri Jan 22 17:19:12 CET 2021
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
>>>
>>> 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.
>>
>>
>> After some testing with a reprex script (included below) we found that
>
> This is broken because email breaks lines unless protected, so this very
> large text example does not run at all. If so much data is needed, and the
> example cannot be retricted to less than 6 points (2D only, no other columns
> needed), post just a script downloading an RDS object (only data frame, no
> tibbles please), convering to the sp object you have problems with, and
> definitely include the spTransform() call (use rgdal directly, report
> rgdal::rgdal_extSoftVersion()).
>
>>
>> 1- No errors or warnings when using the sf package with any recent proj
>> version
>>
>
> Irrelevant, rgdal warns because users need to see that the CRS representation
> has changed. Remember that Proj4 strings should not be used at all, or if
> used will typically be off by 125m (because there is no +datum= or +towgs84=
> mechanism, they are deprecated). The fact that sf doesn't warn is because it
> wasn't thought necessary to alert users to possible output changes. It
> doesn't mean that output doesn't change.
>
>> 2- Switching to the sp package, we found no errors or warnings when the
>> proj version was <=6.3
>>
>> 3- However, the script that uses the sp package, with proj 7.2 does show
>> warnings:
>>
>>
>>> source("aeqd_reprex.R")
>
> Recreating the file by copying from the email to file gives:
>
>> source("aeqd_reprex.R")
> Error in source("aeqd_reprex.R") : aeqd_reprex.R:14:2: unexpected input
> 13: 39.172, 39.172,
> 14:
>
> I'm not going to waste time removing newlines inserted by an email client,
> did you try to extract the fish from the fish soup?
>
> Roger
>
>> [1] "28 Towers; CRS: +proj=longlat +ellps=WGS84 +no_defs"
>> [1] "After transform CRS: +proj=aeqd +lat_0=53.123835 +lon_0=7.03603
>> +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
>> There were 17 warnings (use warnings() to see them)
>>> warnings()[1]
>> Warning message:
>> In showSRID(uprojargs, format = "PROJ", multiline = "NO", ... :
>> Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
>>> sessionInfo()
>> R version 4.0.3 (2020-10-10)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>> Running under: Debian GNU/Linux bullseye/sid
>>
>> Matrix products: default
>> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
>> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so
>>
>> locale:
>> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
>> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
>> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
>> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
>> [9] LC_ADDRESS=C LC_TELEPHONE=C
>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> other attached packages:
>> [1] sp_1.4-4
>>
>> loaded via a namespace (and not attached):
>> [1] compiler_4.0.3 rgdal_1.5-18 grid_4.0.3 lattice_0.20-41
>>
>>
>>> sf::sf_extSoftVersion()[1:3]
>> GEOS GDAL proj.4
>> "3.9.0" "3.2.1" "7.2.0"
>>
>>>
>>> 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
>>>
>>>>
>>>
>>> With a reproducible example that works with proj command line tools, one
>>> could ask for views on the PROJ list itself.
>>>
>>
>>
>> Here is our test script:
>>
>> ##=====================================
>>
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
--
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
More information about the R-sig-Geo
mailing list