Estimating and visualizing python triangulations
Note this is just a proof of concept. We can easily adapt this to multiple animals, make it more streamlined and useful, and tailor it to work easily with whatever format of data you're collecting in the field.
This type of data processing is available in the telemetr package which implements the triangulation approaches detailed in Lenth (1981).
We do not produce error polygons of location estimates, as they're no longer considered very useful. However, the estimation approaches illustrated below do produce standard errors of the X-Y location estimates such that we can visualize uncertainty fairly easily (and could even produce error ellipses with a little more effort).
First, let's make sure we have the necessary R packages and install them if not.
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools", quiet = TRUE)
if (!requireNamespace("telemetr", quietly = TRUE)) devtools::install_github("barryrowlingson/telemetr")
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman", quiet = TRUE)
pacman::p_load(telemetr, dplyr, sf, leaflet)Here's the data you provided:
df <- data.frame(
Date = '10/23/2019',
Observers = 'JoshS',
GID = 1,
Time = c(1,2,3,4),
Easting = c(564821,564830,564839,564846),
Northing = c(2786330,2786317,2786315,2786319),
Azimuth = c(148,60,1,310)
)With telemetr, we have to specify the columns representing the projected X-Y coordinates. We create this as a separate object so we can use the data.frame later.
df_sp <- df
coordinates(df_sp) <- ~ Easting + NorthingNow crudely plot the receiving locations and the observed directions of the signal. One of them gets a little wayward, so we should try a couple of different estimators for the transmitter location, including one that's more robust to the occasional wayward bearing:
plot(df_sp)
telemetr:::drawVectors(~ Azimuth, df_sp)
There are several estimators available to locate the source of the signal. Here we illustrate a maximum likelihood (MLE) approach, which can be susceptible to outliers in observer bearings, and a repeated median regression (RMR) approach that is less sensitive to those outliers.
The estimates are fairly similar, and both provide standard errors of the locational estimate.
# Locate with maximum likelihood
loc_mle <- triang(~ Azimuth, df_sp, method = "mle")
# Locate with repeated median regression
loc_rmr <- triang(~ Azimuth, df_sp, method = "rmr")
# The estimates are similar, and standard errors of the estimate are provided
loc_mle## coordinates ID sd kappa cor se.x se.y ijob err
## 1 (564839.3, 2786322) 1 2.5 525.9321 -0.2222427 0.2753276 0.3176595 1 0
## npts
## 1 4loc_rmr## coordinates ID se.x se.y cor npts
## 1 (564840, 2786324) 1 1.161564 1.428357 -0.07802799 4Plot the resulting estimated locations (blue square == MLE; green triangle == RMR):
# Plot the resulting estimated location
plot(NA, xlim = range(coordinates(df_sp)[,1], coordinates(loc_mle)[,1], coordinates(loc_rmr)[,1]),
ylim = range(coordinates(df_sp)[,2], coordinates(loc_mle)[,2], coordinates(loc_rmr)[,2]),
xlab = "Easting", ylab = "Northing")
points(df_sp)
telemetr:::drawVectors(~ Azimuth, df_sp)
points(loc_mle,pch=0,col="blue")
points(loc_rmr,pch=2,col="green")
These estimates have more meaning when there is some better context provided by, e.g., aerial imagery.
Here are the observation and location estimate data in a quick and dirty interactive map form. In this example, we simply use the maximum directional standard error to illustrate the uncertainty in the location estimates:
# Create a spatial (sf) object of observation locations, assuming WGS84, UTM Zone 17N
df_sf <- st_as_sf(df, coords = c("Easting", "Northing"), crs = 32617)
# Create a spatial (sf) object of observation locations, assuming WGS84, UTM Zone 17N
results <- bind_rows(
as.data.frame(loc_mle),
as.data.frame(loc_rmr)
) %>%
mutate(estimator = c("MLE", "RMR")) %>%
st_as_sf(coords = c("x", "y"), crs = 32617)
# Transform them to geographic coordinates for visualization
df_sf_ll <- st_transform(df_sf, crs = 4326)
results_ll <- st_transform(results, crs = 4326)
# Here we use the maximum SE to illustrate uncertainty in the location estimate
estPal <- colorFactor(palette = c("#008837", "#7b3294"), domain = c("MLE", "RMR"))
p <- leaflet() %>%
addProviderTiles("Esri.WorldImagery",
options = tileOptions(maxNativeZoom=19, maxZoom=21)) %>%
addMarkers(data = df_sf_ll, group = "Telemetry locations") %>%
addCircles(data = results_ll,
radius = ~ se.y, opacity = 0.75, color = ~estPal(estimator),
fillOpacity = 0,
label = ~ estimator) %>%
addLegend("topright", pal = estPal, values = results_ll$estimator,
title = "Estimator", opacity = 0.75) %>%
addScaleBar("bottomright")
pMultiple animal/triangulation example
df2 <- data.frame(
tag = c(rep("A", 4), rep("B", 4), rep("C", 4)),
x = rep(c(281614, 281496, 281439, 281488), 3),
y = rep(c(3753165, 3753217, 3753181, 3753134), 3),
bearing = c(270,175,95,20,265,225,225,275,225,130,110,85),
stringsAsFactors = FALSE
)
actual <- data.frame(
tag = c("A", "B", "C"),
x = c(281507.8, 281402, 281597),
y = c(3753169.4, 3753147.9, 3753152.7),
stringsAsFactors = FALSE
)df2_sp <- df2
coordinates(df2_sp) <- ~ x + y
plot(df2_sp)
telemetr:::drawVectors(~ bearing|tag, df2_sp)
# Locate with maximum likelihood
loc_mle2 <- triang(~ bearing|tag, df2_sp, method = "mle")
# Locate with repeated median regression
loc_rmr2 <- triang(~ bearing|tag, df2_sp, method = "rmr")# Create a spatial (sf) object of observation locations, assuming WGS84, UTM Zone 17N
df2_sf <- st_as_sf(df2, coords = c("x", "y"), crs = 32617)
# Create a spatial (sf) object of observation locations, assuming WGS84, UTM Zone 17N
results2 <- bind_rows(
as.data.frame(loc_mle2),
as.data.frame(loc_rmr2)
) %>%
mutate(estimator = ifelse(!is.na(sd), "MLE", "RMR")) %>%
st_as_sf(coords = c("x", "y"), crs = 32617)
# Write a function to calculate error polygons
error_polys <- function(triang_sf, level = 0.95) {
out <- lapply(seq_len(nrow(triang_sf)), function(i) {
tmp <- triang_sf[i, ]
tmp_coords <- st_coordinates(tmp)
tmp_ellipse <- ellipse::ellipse(tmp$cor, scale = c(tmp$se.x, tmp$se.y),
level = level)
tmp_ellipse[, "x"] <- tmp_ellipse[, "x"] + tmp_coords[, "X"]
tmp_ellipse[, "y"] <- tmp_ellipse[, "y"] + tmp_coords[, "Y"]
tmp_poly <- st_cast(st_multipoint(tmp_ellipse), "POLYGON")
tmp_poly
})
out <- st_as_sfc(out)
out_df <- as.data.frame(triang_sf) %>%
select(tag, cor:se.y, npts, estimator) %>%
mutate(level = level)
out <- st_sf(out_df, out) %>% st_set_crs(st_crs(triang_sf))
out$area <- st_area(out)
out
}
# Calculate error polygons
res_polys <- error_polys(results2)# Transform them to geographic coordinates for visualization
df2_sf_ll <- st_transform(df2_sf, crs = 4326)
results2_ll <- st_transform(results2, crs = 4326)
res_polys_ll <- st_transform(res_polys, crs = 4326)
actual_ll <- st_as_sf(actual, coords = c("x", "y"), crs = 32617) %>%
st_transform(crs = 4326)
# Here we illustrate uncertainty in the location estimate using the 95% error polygons
estPal <- colorFactor(palette = c("#008837", "#7b3294"), domain = c("MLE", "RMR"))
p <- leaflet() %>%
addProviderTiles("Esri.WorldImagery",
options = tileOptions(maxNativeZoom=19, maxZoom=21)) %>%
addAwesomeMarkers(data = filter(df2_sf_ll, tag == "A"),
icon = awesomeIcons("child", library = "fa"),
group = "Telemetry locations")
tags <- unique(df2$tag)
for (t in tags) {
p <- p %>%
addPolygons(data = filter(res_polys_ll, tag == t),
opacity = 0.75, color = ~estPal(estimator),
popup = ~paste(paste("Tag: ", tag),
paste("Estimator:", estimator),
paste("Area:", round(area)),
sep = "<br/>"),
fillOpacity = 0, group = t) %>%
addCircles(data = filter(results2_ll, tag == t),
radius = 1, opacity = 0.75, color = ~estPal(estimator),
fillOpacity = 0, group = t) %>%
addMarkers(data = filter(actual_ll, tag == t), group = t)
}
p <- p %>%
addLegend("topright", pal = estPal, values = unique(results2_ll$estimator),
title = "Estimator", opacity = 0.75) %>%
addScaleBar("bottomright") %>%
leaflet::addLayersControl(overlayGroups = tags,
position = "topleft",
options = leaflet::layersControlOptions(collapsed = FALSE))
p