How to define a local equidistant CRS
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