In this vignette we will analyse the example data from the
moveWindSpeed package. In order to do so we will first load
the data and explore it before diving in to the analysis.
require(moveWindSpeed)## Loading required package: moveWindSpeed## Loading required package: move## Loading required package: geosphere## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0## Loading required package: sp## Loading required package: rasterdata("storks")
plot(
storks,
xlab = "Longitude",
ylab = "Latitude",
main = "Stork data" ,
pch = 19,
cex = .3
)
We focus on 2 minutes where we know that high resolution trajectories were recorded while the bird was in a thermal and select the first individual. We transform the trajectory to a local projection for better plotting.
storksSub <-
storks[strftime(timestamps(storks), format = "%H%M", tz="UTC") %in% c("1303", "1304"), ]
storksWind <- getWindEstimates(storksSub)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /home/bart/.local/share/proj:/usr/share/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.# Use evolution status 2 to avoid using rgdal (set using sp)
set_evolution_status(2L)## [1] 2storksWind <- spTransform(storksWind, center = T)
indiv1 <- storksWind[[1]]
indiv1## class : Move
## features : 120
## extent : -246.2953, 313.526, -55.89234, 45.84754 (xmin, xmax, ymin, ymax)
## crs : +proj=aeqd +lat_0=47.16786075 +lon_0=8.03556635 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## variables : 21
## names : eobs_horizontal_accuracy_estimate, eobs_speed_accuracy_estimate, eobs_status, eobs_temperature, eobs_type_of_fix, eobs_used_time_to_get_fix, ground_speed, heading, height_above_ellipsoid, manually_marked_outlier, event_id, estimationSuccessful, residualVarAirspeed, windX, windY, ...
## min values : 1.54, 0.13, A, 19, 3, 179, 3.63, 1.97, 792.2, NA, 384663102, 1, 0.0849553440610402, 4.89632652004271, -0.0294280996271222, ...
## max values : 1.79, 0.26, A, 19, 3, 299, 16.56, 355.86, 983.7, NA, 384663221, 1, 0.348066492589884, 6.14623034548685, 0.629704086920015, ...
## timestamps : 2014-08-18 13:03:00 ... 2014-08-18 13:04:59 Time difference of 2 mins (start ... end, duration)
## sensors : GPS
## indiv. data : individual_id, tag_id, id, comments, death_comments, earliest_date_born, exact_date_of_birth, latest_date_born, local_identifier, ring_id, sex, taxon_canonical_name
## indiv. value: 21977647 21232176 21977700 4 nestlings. AN873 is the youngest. died 2015-01-11 on landfill Dos Hermanas / S-Spain. Weak strange movements for about 6 hrs. Carcass found, no external injuries. Poisoned? lat: 37.2198845; long: -5.8800427
## NA NA 2014-04-15 00:00:00.000 Bubbel + / DER AN873 (eobs 3666) DER AN873 f NA
## unused rec. : 303
## date created: 2023-06-05 22:37:38Lets first assure that the individual is sampled with one hertz. This is crucial because some of the calculations we do below for plotting assume this.
unique(as.numeric(diff(timestamps(indiv1)), units='secs'))## [1] 1Now lets have a look at the track in one individual within the thermal. We add an arrow to each 10th location to indicate the estimated wind speed.
plot(
indiv1,
type = 'l',
xlim = c(-150, 150),
main = row.names(idData(indiv1))
)
s <- as.numeric(timestamps(indiv1)) %% 10 == 0
expansion <- 3
arrows(
coordinates(indiv1)[s, 1],
coordinates(indiv1)[s, 2],
coordinates(indiv1)[s, 1] + indiv1$windX[s] * expansion,
coordinates(indiv1)[s, 2] + indiv1$windY[s] * expansion,
col = 'red',
length = .1
)
We can also select one rotation and plot how head wind speed vector combined with the airspeed vector results in the ground speed.
indivSub <- indiv1[45:72,]
plot(indivSub, pch = 19)
arrows(
coordinates(indivSub)[-n.locs(indivSub), 1],
coordinates(indivSub)[-n.locs(indivSub), 2],
coordinates(indivSub)[-1, 1],
coordinates(indivSub)[-1, 2],
length = .1
)
arrows(
coordinates(indivSub)[-n.locs(indivSub), 1],
coordinates(indivSub)[-n.locs(indivSub), 2],
coordinates(indivSub)[-n.locs(indivSub), 1] + head(indivSub$windX,-1),
coordinates(indivSub)[-n.locs(indivSub), 2] + head(indivSub$windY,-1),
col = "red",
length = .1
)
arrows(
coordinates(indivSub)[-n.locs(indivSub), 1] + head(indivSub$windX,-1),
coordinates(indivSub)[-n.locs(indivSub), 2] + head(indivSub$windY,-1),
coordinates(indivSub)[-1, 1],
coordinates(indivSub)[-1, 2],
col = "green",
length = .1
)
legend(
"topright",
legend = c("Ground speed", "Wind Speed", "Air Speed"),
bty = "n",
col = c("black", "red", "green"),
lty = 1
)
Now lets see how the track looks once we subtract the wind speed, so that we only look at the rotations relative to the air. We see that the bird half way the thermal changes its relative position in the air column.
s <- !is.na(indiv1$windX)
x <- cumsum(diff(coordinates(indiv1)[, 1])[s] - indiv1$windX[s])
y <- cumsum(diff(coordinates(indiv1)[, 2])[s] - indiv1$windY[s])
plot(x, y, type = 'l', main = "Wind corrected trajectory")
Lets take the same two minute section and compare all wind speed estimates. We see that there all estimates fall in a pretty narrow range.
storksSub <- getWindEstimates(storksSub, windowSize = 31)
plot(
storksSub$windX,
storksSub$windY,
col = trackId(storksSub),
pch = 19,
xlim = range(na.omit(c(0, storksSub$windX))),
ylim = range(na.omit(c(0, storksSub$windY))),
asp = 1,
main = "Distribution of wind speed estimates"
)
We can also create a vertical wind speed profile through a thermal. To do so we select a thermal where most birds are present. We see that there seems to be quite a consistent pattern.
windData <- getWindEstimates(storks)
windData <- windData[strftime(timestamps(windData), format = "%H", tz="UTC") == "12" &
!is.na(windData$windX), ]
plot(
sqrt(windData$windX ^ 2 + windData$windY ^ 2),
windData$height_above_ellipsoid,
main = "Vertical wind profile",
xlab = "Wind speed",
ylab = "Altitude",
col = trackId(windData),
pch = 19
)
We see that if we raise our criteria for thermal soaring the conditions are more rarely occurring.
indivSub <- storks[[1]]
quater <-
getWindEstimates(indivSub, isThermallingFunction = getDefaultIsThermallingFunction(90))
half <-
getWindEstimates(indivSub, isThermallingFunction = getDefaultIsThermallingFunction(180))
full <-
getWindEstimates(indivSub, isThermallingFunction = getDefaultIsThermallingFunction(360))
two <-
getWindEstimates(indivSub, isThermallingFunction = getDefaultIsThermallingFunction(720))
sum(!is.na(quater$windX))## [1] 2551sum(!is.na(half$windX))## [1] 1959sum(!is.na(full$windX))## [1] 803sum(!is.na(two$windX))## [1] 175We can for example also focus on a longer window, meaning we find are more likely to find hits for our rotation criteria.
short <-
getWindEstimates(indivSub,
isThermallingFunction = getDefaultIsThermallingFunction(720),
windowSize = 29)
long <-
getWindEstimates(indivSub,
isThermallingFunction = getDefaultIsThermallingFunction(720),
windowSize = 41)
sum(!is.na(short$windX))## [1] 175sum(!is.na(long$windX))## [1] 338We can also speed up calculations by only looking at each 5th location. This will also reduce the overlap between windows.
system.time(getWindEstimates(
indivSub,
isFocalPoint = function(i, ts) {
T
}
))## user system elapsed
## 4.832 0.696 2.909system.time(getWindEstimates(
indivSub,
isFocalPoint = function(i, ts) {
(i %% 5) == 0
}
))## user system elapsed
## 0.856 0.240 0.574