[R-sig-Geo] How to define a local equidistant CRS
Micha Silver
t@v|b@r @end|ng |rom gm@||@com
Fri Jan 22 08:56:12 CET 2021
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
1- No errors or warnings when using the sf package with any recent proj
version
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")
[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
>>
>>
>
--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
More information about the R-sig-Geo
mailing list