[R-sig-Geo] How to define a local equidistant CRS
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Fri Jan 22 15:30:19 CET 2021
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:
>
> ##=====================================
>
> # library(sf)
> # Option using sp
> library(sp)
> # print(sessionInfo())
> # print(sf_extSoftVersion())
>
> # ---------------------------------------
> # Small subset of microwave data
> # ---------------------------------------
> mw_df = structure(list(Frequency = c(39.172, 37.912, 39.239, 26.362, 39.172,
> 39.172,
> 38.115, 39.165, 37.905, 39.256, 18.03,
> 18.003,
> 38.003, 39.263, 39.263, 38.003, 39.172,
> 37.912,
> 39.239, 39.172, 26.362, 39.172, 38.115,
> 39.165,
> 18.03, 39.256, 37.905, 18.003, 39.165,
> 39.256,
> 37.905, 18.003, 18.03, 38.003, 39.263,
> 39.172,
> 37.912, 39.239, 26.362, 39.172, 39.172,
> 38.115,
> 38.003, 39.263, 39.172, 37.912, 39.239,
> 26.362,
> 39.172, 39.172, 38.115, 37.905, 18.003,
> 39.256,
> 39.165, 18.03, 39.165, 18.03, 37.905,
> 18.003,
> 39.256, 38.003, 39.263, 39.172, 37.912,
> 39.239,
> 26.362, 39.172, 39.172, 38.115, 38.003,
> 39.263,
> 39.165, 39.256, 37.905, 18.03, 18.003,
> 39.172,
> 37.912, 39.239, 26.362, 39.172, 39.172,
> 38.115,
> 38.003, 39.263, 39.172, 39.165, 18.03,
> 37.905,
> 18.003, 39.256, 37.912, 39.239, 26.362,
> 39.172,
> 39.172, 38.115),
> DateTime = c(201209012030, 201209012030, 201209012030,
> 201209012030, 201209012030,
> 201209012030, 201209012030, 201209012030,
> 201209012030, 201209012030,
> 201209012030, 201209012030, 201209012030,
> 201209012030, 201209012045,
> 201209012045, 201209012045, 201209012045,
> 201209012045, 201209012045,
> 201209012045, 201209012045, 201209012045,
> 201209012045, 201209012045,
> 201209012045, 201209012045, 201209012045,
> 201209012100, 201209012100,
> 201209012100, 201209012100, 201209012100,
> 201209012100, 201209012100,
> 201209012100, 201209012100, 201209012100,
> 201209012100, 201209012100,
> 201209012100, 201209012100, 201209012115,
> 201209012115, 201209012115,
> 201209012115, 201209012115, 201209012115,
> 201209012115, 201209012115,
> 201209012115, 201209012115, 201209012115,
> 201209012115, 201209012115,
> 201209012115, 201209012130, 201209012130,
> 201209012130, 201209012130,
> 201209012130, 201209012130, 201209012130,
> 201209012130, 201209012130,
> 201209012130, 201209012130, 201209012130,
> 201209012130, 201209012130,
> 201209012145, 201209012145, 201209012145,
> 201209012145, 201209012145,
> 201209012145, 201209012145, 201209012145,
> 201209012145, 201209012145,
> 201209012145, 201209012145, 201209012145,
> 201209012145, 201209012200,
> 201209012200, 201209012200, 201209012200,
> 201209012200, 201209012200,
> 201209012200, 201209012200, 201209012200,
> 201209012200, 201209012200,
> 201209012200, 201209012200, 201209012200),
> Pmin = c(-59, -52, -55, -49, -50, -51, -45, -46, -48,
> -54, -59, -56, -45, -44, -44,
> -45, -59, -52, -55, -50, -49, -51, -45, -46,
> -59, -54, -48, -56, -47, -54,
> -48, -56, -58, -44, -44, -59, -52, -56, -49,
> -50, -51, -45, -44, -44, -60,
> -52, -56, -50, -51, -50, -45, -48, -56, -54,
> -47, -57, -46, -57, -48, -56,
> -54, -44, -44, -60, -52, -56, -50, -51, -50,
> -45, -44, -44, -47, -54, -48,
> -58, -57, -60, -52, -55, -49, -50, -51, -45,
> -44, -45, -60, -45, -58, -48,
> -57, -54, -52, -55, -50, -51, -50, -45),
> Pmax = c(-58, -52, -54, -46, -49, -51, -45, -44, -47,
> -53, -57, -55, -44, -43, -43,
> -44, -58, -52, -54, -49, -46, -51, -45, -44,
> -57, -53, -47, -55, -46, -53,
> -47, -55, -57, -43, -43, -58, -52, -54, -46,
> -49, -51, -45, -44, -43, -58,
> -52, -54, -47, -51, -49, -44, -47, -55, -53,
> -43, -57, -43, -57, -47, -55,
> -53, -44, -43, -58, -52, -54, -47, -51, -49,
> -45, -44, -44, -44, -53, -47,
> -57, -55, -59, -52, -54, -46, -49, -51, -45,
> -44, -44, -59, -44, -57, -47,
> -56, -53, -52, -54, -47, -51, -49, -44),
> PathLength = c(2.97052, 1.50853, 5.97957, 8.19688,
> 2.26532, 1.50853, 1.504, 4.78188,
> 4.78188, 2.12735, 13.50599, 7.56082, 4.37035,
> 4.37035, 4.37035,
> 4.37035, 2.97052, 1.50853, 5.97957, 2.26532,
> 8.19688, 1.50853,
> 1.504, 4.78188, 13.50599, 2.12735, 4.78188,
> 7.56082, 4.78188, 2.12735,
> 4.78188, 7.56082, 13.50599, 4.37035, 4.37035,
> 2.97052, 1.50853,
> 5.97957, 8.19688, 2.26532, 1.50853, 1.504,
> 4.37035, 4.37035, 2.97052,
> 1.50853, 5.97957, 8.19688, 1.50853, 2.26532,
> 1.504, 4.78188, 7.56082,
> 2.12735, 4.78188, 13.50599, 4.78188, 13.50599,
> 4.78188, 7.56082,
> 2.12735, 4.37035, 4.37035, 2.97052, 1.50853,
> 5.97957, 8.19688,
> 1.50853, 2.26532, 1.504, 4.37035, 4.37035,
> 4.78188, 2.12735, 4.78188,
> 13.50599, 7.56082, 2.97052,1.50853, 5.97957,
> 8.19688, 2.26532,
> 1.50853, 1.504, 4.37035, 4.37035, 2.97052,
> 4.78188, 13.50599, 4.78188,
> 7.56082, 2.12735, 1.50853, 5.97957, 8.19688,
> 1.50853, 2.26532, 1.504),
> XStart = c(6.98795,6.91965, 7.12605, 7.12605, 6.90439,
> 6.92231, 7.05273, 6.94839,
> 7.01525, 6.95347, 6.94267, 7.12769, 6.92516,
> 6.98553, 6.98553, 6.92516,
> 6.98795, 6.91965, 7.12605, 6.90439, 7.12605,
> 6.92231, 7.05273, 6.94839,
> 6.94267, 6.95347, 7.01525, 7.12769, 6.94839,
> 6.95347, 7.01525, 7.12769,
> 6.94267, 6.92516, 6.98553, 6.98795, 6.91965,
> 7.12605, 7.12605, 6.90439,
> 6.92231, 7.05273, 6.92516, 6.98553, 6.98795,
> 6.91965, 7.12605, 7.12605,
> 6.92231, 6.90439,7.05273, 7.01525, 7.12769,
> 6.95347, 6.94839, 6.94267,
> 6.94839, 6.94267, 7.01525, 7.12769, 6.95347,
> 6.92516, 6.98553, 6.98795,
> 6.91965, 7.12605, 7.12605, 6.92231, 6.90439,
> 7.05273, 6.92516, 6.98553,
> 6.94839, 6.95347, 7.01525, 6.94267, 7.12769,
> 6.98795, 6.91965, 7.12605,
> 7.12605, 6.90439, 6.92231, 7.05273, 6.92516,
> 6.98553, 6.98795, 6.94839,
> 6.94267, 7.01525, 7.12769, 6.95347, 6.91965,
> 7.12605, 7.12605, 6.92231,
> 6.90439, 7.05273),
> YStart = c(53.09842, 53.31928, 53.14836, 53.14836,
> 53.32724, 53.33274, 53.1538,
> 52.91493, 52.92956, 52.98836, 52.99448, 52.92908,
> 53.17983, 53.16479,
> 53.16479, 53.17983, 53.09842, 53.31928, 53.14836,
> 53.32724, 53.14836,
> 53.33274, 53.1538, 52.91493, 52.99448, 52.98836,
> 52.92956, 52.92908,
> 52.91493, 52.98836, 52.92956, 52.92908, 52.99448,
> 53.17983, 53.16479,
> 53.09842, 53.31928, 53.14836, 53.14836, 53.32724,
> 53.33274, 53.1538,
> 53.17983, 53.16479, 53.09842, 53.31928, 53.14836,
> 53.14836, 53.33274,
> 53.32724, 53.1538, 52.92956, 52.92908, 52.98836,
> 52.91493, 52.99448,
> 52.91493, 52.99448, 52.92956, 52.92908, 52.98836,
> 53.17983, 53.16479,
> 53.09842, 53.31928, 53.14836, 53.14836, 53.33274,
> 53.32724, 53.1538,
> 53.17983, 53.16479, 52.91493, 52.98836, 52.92956,
> 52.99448, 52.92908,
> 53.09842, 53.31928, 53.14836, 53.14836, 53.32724,
> 53.33274, 53.1538,
> 53.17983, 53.16479, 53.09842, 52.91493, 52.99448,
> 52.92956, 52.92908,
> 52.98836, 53.31928, 53.14836, 53.14836, 53.33274,
> 53.32724, 53.1538),
> XEnd = c(6.948, 6.92231, 7.20104, 7.02248, 6.87102, 6.91965,
> 7.05556, 7.01525,
> 6.94839, 6.97979, 7.07683, 7.01525, 6.98553, 6.92516,
> 6.92516, 6.98553,
> 6.948, 6.92231, 7.20104, 6.87102, 7.02248, 6.91965,
> 7.05556, 7.01525,
> 7.07683, 6.97979, 6.94839, 7.01525, 7.01525, 6.97979,
> 6.94839, 7.01525,
> 7.07683, 6.98553, 6.92516, 6.948, 6.92231, 7.20104,
> 7.02248, 6.87102,
> 6.91965, 7.05556, 6.98553, 6.92516, 6.948, 6.92231,
> 7.20104, 7.02248,
> 6.91965, 6.87102, 7.05556, 6.94839, 7.01525, 6.97979,
> 7.01525, 7.07683,
> 7.01525, 7.07683, 6.94839, 7.01525, 6.97979, 6.98553,
> 6.92516, 6.948,
> 6.92231, 7.20104, 7.02248, 6.91965, 6.87102, 7.05556,
> 6.98553, 6.92516,
> 7.01525, 6.97979, 6.94839, 7.07683, 7.01525, 6.948,
> 6.92231, 7.20104,
> 7.02248, 6.87102, 6.91965, 7.05556, 6.98553, 6.92516,
> 6.948, 7.01525,
> 7.07683, 6.94839, 7.01525, 6.97979, 6.92231, 7.20104,
> 7.02248, 6.91965,
> 6.87102, 7.05556),
> YEnd = c(53.08685, 53.33274, 53.17761, 53.10906,
> 53.32334, 53.31928, 53.14039,
> 52.92956, 52.91493, 52.97772, 53.08498, 52.92956,
> 53.16479, 53.17983,
> 53.17983, 53.16479, 53.08685, 53.33274, 53.17761,
> 53.32334, 53.10906,
> 53.31928, 53.14039, 52.92956, 53.08498, 52.97772,
> 52.91493, 52.92956,
> 52.92956, 52.97772, 52.91493, 52.92956, 53.08498,
> 53.16479, 53.17983,
> 53.08685, 53.33274, 53.17761, 53.10906, 53.32334,
> 53.31928, 53.14039,
> 53.16479, 53.17983, 53.08685, 53.33274, 53.17761,
> 53.10906, 53.31928,
> 53.32334, 53.14039, 52.91493, 52.92956, 52.97772,
> 52.92956, 53.08498,
> 52.92956, 53.08498, 52.91493, 52.92956, 52.97772,
> 53.16479, 53.17983,
> 53.08685, 53.33274, 53.17761, 53.10906, 53.31928,
> 53.32334, 53.14039,
> 53.16479, 53.17983, 52.92956, 52.97772, 52.91493,
> 53.08498, 52.92956,
> 53.08685, 53.33274, 53.17761, 53.10906, 53.32334,
> 53.31928, 53.14039,
> 53.16479, 53.17983, 53.08685, 52.92956, 53.08498,
> 52.91493, 52.92956,
> 52.97772, 53.33274, 53.17761, 53.10906, 53.31928,
> 53.32334, 53.14039),
> ID = c(901, 926, 935, 939, 936, 937, 1030, 471, 543, 478, 475,
> 546, 818, 820,
> 820, 818, 901, 926, 935, 936, 939, 937, 1030, 471, 475,
> 478, 543, 546,
> 471, 478, 543, 546, 475, 818, 820, 901, 926, 935, 939,
> 936, 937, 1030,
> 818, 820, 901, 926, 935, 939, 937, 936, 1030, 543, 546,
> 478, 471, 475,
> 471, 475, 543, 546, 478, 818, 820, 901, 926, 935, 939,
> 937, 936, 1030,
> 818, 820, 471, 478, 543, 475, 546, 901, 926, 935, 939,
> 936, 937, 1030,
> 818, 820, 901, 471, 475, 543, 546, 478, 926, 935, 939,
> 937, 936, 1030)),
> class = c("spec_tbl_df", "tbl_df", "tbl", "data.frame"),
> row.names = c(NA, -98L),
> spec = structure(list(cols = list(Frequency = structure(list(),
> class = c("collector_double", "collector")),
> DateTime = structure(list(),
> class = c("collector_double", "collector")),
> Pmin = structure(list(), class =
> c("collector_double", "collector")),
> Pmax = structure(list(), class =
> c("collector_double", "collector")),
> PathLength = structure(list(),
> class = c("collector_double", "collector")),
> XStart = structure(list(), class
> = c("collector_double", "collector")),
> YStart = structure(list(), class
> = c("collector_double", "collector")),
> XEnd = structure(list(), class =
> c("collector_double", "collector")),
> YEnd = structure(list(), class =
> c("collector_double", "collector")),
> ID = structure(list(), class =
> c("collector_double", "collector"))),
> default = structure(list(), class =
> c("collector_guess", "collector")),
> skip = 1L), class = "col_spec"))
>
> # ---------------------------------------
> # Find middle point
> # Setup CRS strings
> # ---------------------------------------
> XMiddle <- (min(mw_df$XStart, mw_df$XEnd) + max(mw_df$XStart, mw_df$XEnd)) /
> 2
> YMiddle <- (min(mw_df$YStart, mw_df$YEnd) + max(mw_df$YStart, mw_df$YEnd)) /
> 2
>
> projstring_wgs84 = "+proj=longlat +ellps=WGS84"
> projstring_aeqd <- paste("+proj=aeqd +lat_0=", YMiddle,
> " +lon_0=", XMiddle,
> " +x_0=0 +y_0=0 +ellps=WGS84",sep="")
>
> # ---------------------------------------
> # Loop thru all ID's and extract both start and end points
> # Save in list, to get spatial points object of all MW towers
> # ---------------------------------------
> IDs = unique(mw_df$ID)
> towers_list = lapply(IDs, function(i) {
> # Create point sf object with two features, start and end of link
> link = mw_df[mw_df$ID == i,]
> # All time slots have the same coords. Get just the first:
> x1 = link$XStart[1]
> y1 = link$YStart[1]
> x2 = link$XEnd[1]
> y2 = link$YEnd[1]
> towers_df = data.frame("x" = c(x1, x2), "y" = c(y1, y2),
> "ID"=i, "End" = c("Start", "End"))
>
> # Create sf object of point features
> #towers_sf = st_as_sf(towers_df, coords = c("x", "y"),
> # crs = projstring_wgs84, agr = "constant")
> #return(towers_sf)
> # Option using sp
> xy = towers_df[,c("x", "y")]
> towers_sp = SpatialPointsDataFrame(coords = xy, data = towers_df,
> proj4string = CRS(projstring_wgs84))
> return(towers_sp)
> })
>
> # ---------------------------------------
> # Bind all point together and transform to AEQD
> # ---------------------------------------
> towers = do.call(rbind, towers_list)
> num_towers = length(towers$ID)
>
> # sf Option
> # ---------------------------------------
> # st_crs(towers) = projstring_wgs84
> # print(paste(num_towers, "Towers; CRS:", st_crs(towers)$proj4string))
> # towers_aeqd = st_transform(towers, projstring_aeqd)
> # print(paste("After transform CRS:", st_crs(towers_aeqd)$proj4string))
>
> # sp Option
> # ---------------------------------------
> print(paste(num_towers, "Towers; CRS:", proj4string(towers)))
> towers_aeqd = spTransform(towers, CRS(projstring_aeqd))
> print(paste("After transform CRS:", proj4string(towers_aeqd)))
>
> ##=====================================
>
>> Hope this helps,
>>
>> Roger
>>
>
> Thanks for your comments,
>
> Micha
>
>
>
>>> 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
More information about the R-sig-Geo
mailing list