R/new_plotting_functions.R
In eldarrak/FLightR: Reconstruct Animal Paths from Solar Geolocation Loggers Data
Defines functions
get_utilisation_points
get_points_distribution
plot_likelihood
seasonal_donut
plot_util_distr
get_time_spent_buffer
get_gunion_r
get_buffer
plot_lon_lat
map.FLightR.ggmap
Documented in
map.FLightR.ggmap
plot_likelihood
plot_lon_lat
plot_util_distr
#' plots result over map
#'
#' plots track over map with probability cloud. Can plot only part of the track if dates are specified. Note that you can use it only after obtaining and registering in you current session Google Api Key. For details on the API key check [here](http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error).
#' @param Result FLightR result object obtained from \code{\link{run.particle.filter}}
#' @param dates either NULL if all twilights should be included or data.frame with first column - start of the period and second end of the period. Each line represents a new period
#' @param plot.cloud Should probability cloud be plotted? If TRUE cloud is estimated by \code{\link[ggplot2]{stat_density2d}}
#' @param map.options options passed to \code{\link[ggmap]{get_map}}, note that \code{zoom} option is defined separately
#' @param plot.options plotting options. Not defined yet!
#' @param save.options options passed to \code{\link[ggplot2]{ggsave}}. Filename should be defined here.
#' @param zoom Zoom for map. If 'auto' FLightR will try to find optimal zoom level by downloading different size maps and checking whether all the points fit the map.
#' @param return.ggobj Should ggobj be returned for subsequent checks and/or replotting
#' @param seasonal.colors if true points of the track will have seasonal colors
#' @param seasonal.donut.location if NULL - no color wheel placed, otherwise select one of 'bottomleft', 'bottomright', 'topleft'
#' @param seasonal.donut.proportion how much of X axis should color wheel occupy.
#' return either NULL or ggplot2 class object
#' @param save should function save results with \code{\link[ggplot2]{ggsave}}?
#' @return if 'return.ggobj=TRUE' return ggplot object otherwise returns 'NULL'.
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-06-25', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#' # here we will run only 1e4 partilces for a very short track.
#' # One should use 1e6 particles for the full run
#' Result<-run.particle.filter(all.in, threads=1,
#' nParticles=1e3, known.last=TRUE,
#' precision.sd=25, check.outliers=FALSE)
#'
#'\dontrun{
#' map.FLightR.ggmap(Result, seasonal.donut.location=NULL, zoom=6, save=FALSE)
#'}
#' # for this short track without variance seasonal donut does not work,
#' # but for normall track it will.
#' @author Eldar Rakhimberdiev
#' @export map.FLightR.ggmap
map.FLightR.ggmap<-function(Result, dates=NULL, plot.cloud=TRUE, map.options=NULL, plot.options=NULL, save.options=NULL, zoom="auto", return.ggobj=FALSE, seasonal.colors=TRUE, seasonal.donut.location='topleft', seasonal.donut.proportion=0.5, save=TRUE) {
if (!is.null(plot.options)) warning("plot options are not in use yet. Let me know what you would like to have here.")
if (utils::packageVersion('ggmap')[1]<2.7) { stop('map.FLightR.ggmap function works only with ggmap >= 2.7.x') }
if (!ggmap::has_google_key()) stop('From August 2018 Google allows to use Google maps only for users with the API key, please get one and proceed as described here: http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error/')
# dates should be a data.frame with first point - starting dates and last column end dates for periods
# ggsave.options is a list that will be will be directly passed to ggsave
Opt<-options('ggmap')
if(is.null(Opt$ggmap$google$account_type)) Opt$ggmap$google$account_type<-'standard'
if(is.null(Opt$ggmap$google$day_limit)) Opt$ggmap$google$day_limit<-2500
if(is.null(Opt$ggmap$google$second_limit)) Opt$ggmap$google$day_limit<-50
if(is.null(Opt$ggmap$google$client)) Opt$ggmap$google$client<-NA
if(is.null(Opt$ggmap$google$signature)) Opt$ggmap$google$signature<-NA
options('ggmap'= Opt$ggmap)
# select twilights to plot
if (is.null(dates)) {
twilights.index<-1:length(Result$Results$Points.rle)
} else {
twilights.index<-c()
for (segment in 1:nrow(dates)) {
twilights.index<-c(twilights.index, which(Result$Results$Quantiles$time>=dates[segment,1] & Result$Results$Quantiles$time<=dates[segment,2]))
}
if (length(twilights.index)==0) stop("dates do not overlap with the track time span!")
}
if (plot.cloud) {
Points_rle<-Result$Results$Points.rle[twilights.index]
All.Points<-rep(0, nrow(Result$Spatial$Grid))
for (twilight in 1:length(twilights.index)) {
All.Points[Points_rle[[twilight]]$values]<-All.Points[Points_rle[[twilight]]$values] + Points_rle[[twilight]]$lengths
}
Points<-Result$Spatial$Grid[All.Points>0,,drop=FALSE][sample.int(length(All.Points[All.Points>0]), size = 10000, replace = TRUE, prob = All.Points[All.Points>0]), 1:2]
}
# background map
# check whether Grid was over dateline:
overdateline<-ifelse(attr(Result$Spatial$Grid, 'left')> attr(Result$Spatial$Grid, 'right'), TRUE, FALSE)
if (overdateline) {
location<-cbind(
min(Result$Results$Quantiles$Medianlon[twilights.index][Result$Results$Quantiles$Medianlon[twilights.index]>0]),
min(Result$Results$Quantiles$Medianlat[twilights.index]),
max(Result$Results$Quantiles$Medianlon[twilights.index][Result$Results$Quantiles$Medianlon[twilights.index]<0]),
max(Result$Results$Quantiles$Medianlat[twilights.index]))
if (!is.finite(location[1])) location[1] <-180
if (!is.finite(location[3])) location[3] <- -180
} else {
location<-cbind(
min(Result$Results$Quantiles$Medianlon[twilights.index]),
min(Result$Results$Quantiles$Medianlat[twilights.index]),
max(Result$Results$Quantiles$Medianlon[twilights.index]),
max(Result$Results$Quantiles$Medianlat[twilights.index]))
}
# there could be problems here if the dateline is passed. I will add check for it later.
center_lon<-mean(c(location[1], location[3]))
if (overdateline) center_lon<-mean(c(location[1], location[3]+360))-360
if (center_lon < (-180)) center_lon<-center_lon+360
center_lat<-mean(c(location[2], location[4]))
if (center_lon>180 ) center_lon<-center_lon-360
if ((location[1])>180 ) location[1]<-location[1]-360
if ((location[3])>180 ) location[3]<-location[3]-360
# combine location with map options
if (is.null(map.options)) map.options<-list()
if (is.null(map.options$location)) map.options$location<-c(center_lon, center_lat)
if (is.null(map.options$zoom) & is.numeric(zoom)) map.options$zoom=zoom
if (is.null(map.options$col)) map.options$col="bw"
if (zoom=="auto") {
for (zoom_cur in (2:10)) {
map.options$zoom=zoom_cur
background <-do.call(ggmap::get_map, map.options)
bb<-attr(background, 'bb')
bb[2]<-ifelse(bb[2]< (-180), bb[2]+360, bb[2])
bb[4]<-ifelse(bb[4]< (-180), bb[4]+360, bb[4])
bb[2]<-ifelse(bb[2]> (180), bb[2]-360, bb[2])
bb[4]<-ifelse(bb[4]> (180), bb[4]-360, bb[4])
if (bb[4]<bb[2]) {
lonisinbb<-(location[1] >= bb[2] || location[1] <= bb[4]) & (location[3] >= bb[2] || location[3] <= bb[4])
} else {
lonisinbb<-(location[1] >= bb[2] && location[1] <= bb[4]) & (location[3] >= bb[2] && location[3] <= bb[4])
}
isinbb<-lonisinbb & location[2]>=bb[1] & location[4]<=bb[3]
if (!isinbb) break
}
map.options$zoom=zoom_cur-1
}
background <-do.call(ggmap::get_map, map.options)
p<-ggmap::ggmap(background)#, maprange=TRUE)
if (plot.cloud) {
#if (overdateline) Points[,1]<-ifelse(Points[,1]>0, Points[,1]-360, Points[,1])
# for r cmd check:
..level..<-NA
lon<-NA
lat<-NA
#################
p<- p+ ggplot2::stat_density2d(data=data.frame(Points), ggplot2::aes(fill = ..level.., alpha = ..level.., x=lon, y=lat), size = 0.01, geom = 'polygon', n=400)
p<-p+ ggplot2::scale_fill_gradient(low = "green", high = "red")
p<-p+ ggplot2::scale_alpha(range = c(0.00, 0.25), guide = FALSE)
}
#p<-p+ggplot2::coord_map(projection="mercator",
# xlim=c(attr(background, "bb")$ll.lon, attr(background, "bb")$ur.lon),
# ylim=c(attr(background, "bb")$ll.lat, attr(background, "bb")$ur.lat))
p<-p+ ggplot2::theme(legend.position = "none", axis.title = ggplot2::element_blank(), text = ggplot2::element_text(size = 12))
# here I plot track for selected dates
Points2plot<-data.frame(lat=Result$Results$Quantiles$Medianlat,lon=Result$Results$Quantiles$Medianlon)
if (overdateline) Points2plot[,2][Points2plot[,2]>attr(background, "bb")$ur.lon] <- Points2plot[,2][Points2plot[,2]>attr(background, "bb")$ur.lon]-360
if (overdateline) Points2plot[,2][Points2plot[,2]<attr(background, "bb")$ll.lon] <- Points2plot[,2][Points2plot[,2]<attr(background, "bb")$ll.lon] + 360
if (is.null(dates)) {
cur_twilights<-1:length(Result$Results$Quantiles$time)
Total_segments<-1
} else {
cur_twilights<-which(Result$Results$Quantiles$time>=dates[segment,1] & Result$Results$Quantiles$time<=dates[segment,2])
Total_segments<-nrow(dates)
}
for (segment in 1:Total_segments) {
p<-p+ ggplot2::geom_path(data=Points2plot[cur_twilights,],ggplot2::aes(x=lon,y=lat), colour=grDevices::grey(0.3))
if (seasonal.colors) {
Seasonal_palette<-grDevices::colorRampPalette(grDevices::hsv(1-((1:365)+(365/4))%%365/365, s=0.8, v=0.8), space="Lab")
Days<-as.integer(format(Result$Results$Quantiles$time[cur_twilights], "%j"))
p<-p+ggplot2::geom_point(data=Points2plot[cur_twilights,], shape=21, colour=grDevices::grey(0.3), bg=Seasonal_palette(365)[Days])
} else {
p<-p+ggplot2::geom_point(data=Points2plot[cur_twilights,], shape="+", colour=grDevices::grey(0.3))
}
}
if (!is.null(seasonal.donut.location)) {
Xrange<-c(attr(background, "bb")$ll.lon, attr(background, "bb")$ur.lon)
Yrange<-c(attr(background, "bb")$ll.lat, attr(background, "bb")$ur.lat)
d<-seasonal_donut()
g = ggplot2::ggplotGrob(d)
if (seasonal.donut.location=='bottomleft') {
p<-p + ggmap::inset(grob = g, xmin = Xrange[1], xmax = Xrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]), ymin = Yrange[1], ymax = Yrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]))
}
if (seasonal.donut.location=='topleft') {
p<-p + ggmap::inset(grob = g, xmin = Xrange[1], xmax = Xrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]), ymin = Yrange[2]-seasonal.donut.proportion*(Yrange[2]-Yrange[1]), ymax = Yrange[2])
}
if (seasonal.donut.location=='bottomright') {
p<-p + ggmap::inset(grob = g, xmin = Xrange[2]-seasonal.donut.proportion*(Yrange[2]-Yrange[1]), xmax = Xrange[2], ymin = Yrange[1], ymax = Yrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]))
}
}
graphics::plot(p)
if (save) {
if (is.null(save.options)) save.options<-list()
if (is.null(save.options$filename)) save.options$filename<-"FLightR.map.pdf"
do.call(ggplot2::ggsave, save.options)
}
if (return.ggobj) {return(p)} else {return(NULL)}
}
#' plots result by longitude and latitude
#'
#' This function plots result by latitude and longitude in either vertical or horizontal layout.
#' @param Result FLightR result object obtained from \code{\link{run.particle.filter}}
#' @param scheme either 'vertical' or 'horizontal' layouts
#' return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#' # here we will run only 1e4 partilces for a very short track.
#' # One should use 1e6 particles for the full run
#' Result<-run.particle.filter(all.in, threads=1,
#' nParticles=1e3, known.last=TRUE,
#' precision.sd=25, check.outliers=FALSE)
#'
#' plot_lon_lat(Result)
#'
#' @author Eldar Rakhimberdiev
#' @export plot_lon_lat
plot_lon_lat<-function(Result, scheme=c("vertical", "horizontal")) {
oldpar <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(oldpar))
Quantiles<-Result$Results$Quantiles
# check whether Grid was over dateline:
overdateline<-ifelse(attr(Result$Spatial$Grid, 'left')> attr(Result$Spatial$Grid, 'right'), TRUE, FALSE)
if (overdateline) {
Quantiles$Medianlon[Quantiles$Medianlon<0]<-Quantiles$Medianlon[Quantiles$Medianlon<0]+360
Quantiles$LCI.lon[Quantiles$LCI.lon<0]<-Quantiles$LCI.lon[Quantiles$LCI.lon<0]+360
Quantiles$UCI.lon[Quantiles$UCI.lon<0]<-Quantiles$UCI.lon[Quantiles$UCI.lon<0]+360
Quantiles$TrdQu.lon[Quantiles$TrdQu.lon<0]<-Quantiles$TrdQu.lon[Quantiles$TrdQu.lon<0]+360
Quantiles$FstQu.lon[Quantiles$FstQu.lon<0]<-Quantiles$FstQu.lon[Quantiles$FstQu.lon<0]+360
}
if (scheme[1]=="horizontal") {
graphics::par(mfrow=c(1,2))
} else {
graphics::par(mfrow=c(2,1))
}
graphics::par(mar=c(2,4,3,1),cex=1)
suppressWarnings(Sys.setlocale("LC_ALL", "English"))
#------here I have create the equinox dates..
Years<-unique(format(Quantiles$time, format="%Y"))
eq<-c(as.POSIXct(paste(Years, "-09-22 00:00:00 GMT", sep="")), as.POSIXct(paste(Years, "-03-20 12:00:00 GMT", sep="")))
eq<-eq[eq>min(Quantiles$time) & eq<max(Quantiles$time)]
#------- vert grid
Vert_grid<-seq(as.POSIXct("2000-01-01", tz='GMT'), as.POSIXct("2030-01-01", tz='GMT'), by="month")
Vert_grid<-Vert_grid[Vert_grid>=min(Quantiles$time) & Vert_grid<=max(Quantiles$time)]
#Longitude
graphics::plot(Quantiles$Medianlon~Quantiles$time, las=1,col=grDevices::grey(0.1),pch=16,
ylab="Longitude",xlab="",lwd=2, ylim=range(c( Quantiles$LCI.lon,
Quantiles$UCI.lon )), type="n", axes=FALSE)
graphics::axis(2, las=1)
graphics::axis.POSIXct(1, x=Quantiles$time, format="1-%b")
graphics::box()
# add vertical lines for the first day of every month
graphics::abline(v=Vert_grid, col=grDevices::grey(0.5), lty=2)
graphics::abline(h=seq(-180, 360, by=10), col=grDevices::grey(0.5), lty=2)
graphics::abline(v=eq, col=2, lwd=2, lty=1)
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)),
y=c(Quantiles$LCI.lon, rev(Quantiles$UCI.lon)),
col=grDevices::grey(0.9), border=grDevices::grey(0.5))
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)),
y=c(Quantiles$TrdQu.lon, rev(Quantiles$FstQu.lon)),
col=grDevices::grey(0.7), border=grDevices::grey(0.5))
graphics::lines(Quantiles$Medianlon~Quantiles$time, col=grDevices::grey(0.1),lwd=2)
#Latitude
graphics::par(mar=c(3,4,1,1))
graphics::plot(Quantiles$Medianlat~Quantiles$time, las=1,col=grDevices::grey(0.1),
pch=16,ylab="Latitude",xlab="",lwd=2,
ylim=range(c( Quantiles$UCI.lat, Quantiles$LCI.lat )), type="n", axes=FALSE)
graphics::axis(2, las=1)
graphics::axis.POSIXct(1, x=Quantiles$time, format="1-%b")
graphics::box()
graphics::abline(v=Vert_grid, col=grDevices::grey(0.5), lty=2)
graphics::abline(h=seq(-80, 80, by=10), col=grDevices::grey(0.5), lty=2)
graphics::abline(v=eq, col=2, lwd=2, lty=1)
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)), y=c(Quantiles$LCI.lat, rev(Quantiles$UCI.lat)),
col=grDevices::grey(0.9), border=grDevices::grey(0.5))
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)), y=c(Quantiles$TrdQu.lat, rev(Quantiles$FstQu.lat)),
col=grDevices::grey(0.7), border=grDevices::grey(0.5))
graphics::lines(Quantiles$Medianlat~Quantiles$time, col=grDevices::grey(0.1),lwd=2)
return(NULL)
}
get_buffer<-function(coords, r){
# p = sp::SpatialPoints(matrix(coords, ncol=2), proj4string=sp::CRS("+proj=longlat +datum=WGS84"))
# aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
# p@coords[[2]], p@coords[[1]])
# projected <- sp::spTransform(p, sp::CRS(aeqd))
# buffered <- rgeos::gBuffer(projected, width=r, byid=TRUE)
# buffer_lonlat <- sp::spTransform(buffered, CRS=p@proj4string)
p <- sf::st_as_sf(as.data.frame(matrix(coords, ncol=2)), coords=c('V1','V2'),crs = 4326)
aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
sf::st_coordinates(p)[[2]], sf::st_coordinates(p)[[1]])
projected <- sf::st_transform(p,sf::st_crs(aeqd))
buffered <- sf::st_buffer(projected, dist = r)
buffer_lonlat <- sf::st_transform(buffered,sf::st_crs(p))
return(buffer_lonlat)
}
get_gunion_r<-function(Result) {
#Distances= sp::spDists(Result$Spatial$Grid[1:min(c(nrow(Result$Spatial$Grid), 1000)),1:2], longlat=TRUE)
Distances= as.numeric(sf::st_distance(sf::st_as_sf(as.data.frame(Result$Spatial$Grid[1:min(c(nrow(Result$Spatial$Grid), 1000)),1:2]),
coords=c('lon','lat'),
crs=4326))/1000)
# ok, distances go up to 51.2.. the next step is 62..
# so if I round them
Selected_dist<-unique(sort(round(Distances/10)*10))[2]
r<-Selected_dist*0.85*1000
return(r)
}
get_time_spent_buffer<-function(Result, dates=NULL, percentile=0.5, r=NULL) {
# r in meters..
# dates could be either NULL - for plotting all available dates,
# or numeric with length of one - to plot a specific twilight by number (row number in Quantiles)
# or data.frame with two columns - one for the first date in period and one for the last date in period.
if (is.null(dates)) {
twilights.index<-1:length(Result$Results$Points.rle)
} else {
if (is.numeric(dates) & length(dates)==1) {
twilights.index<-dates
} else {
twilights.index<-c()
for (segment in 1:nrow(dates)) {
twilights.index<-c(twilights.index, which(Result$Results$Quantiles$time>=dates[segment,1] & Result$Results$Quantiles$time<=dates[segment,2]))
}
if (length(twilights.index)==0) stop("dates do not overlap with the track time span!")
}
}
message('function will plot ', length(twilights.index), ' twilights\n')
Points_selected<-get_utilisation_points(Result, twilights.index, percentile)
if (is.null(r)) {
r=get_gunion_r(Result)
}
#Now I want to create a spatial buffers around points with this radius
#spoints = SpatialPoints(Result$Spatial$Grid[Points_selected,1:2], proj4string= CRS("+proj=longlat +datum=WGS84"))
Buffers<-apply(matrix(Result$Spatial$Grid[Points_selected,1:2], ncol=2), 1, get_buffer, r=r)
Buff_comb<-Buffers[[1]]
if (length(Buffers)>1) {
for (i in 2:length(Buffers)) {
#Buff_comb<-rgeos::gUnion(Buff_comb, Buffers[[i]])
Buff_comb<-sf::st_union(Buff_comb, Buffers[[i]])
}
}
#Buff_comb_simpl<-rgeos::gSimplify(Buff_comb, tol=0.01, topologyPreserve=TRUE)
Buff_comb_simpl<-sf::st_simplify(Buff_comb,preserveTopology = FALSE, dTolerance = 0.01)
return(list(Buffer=Buff_comb_simpl, nPoints=length(Points_selected))) # Points before
}
#' plots resulting track over map with uncertainty shown by space utilisation distribution
#'
#' May be use not only for the whole track but for a set of specific dates, e.g. to show spatial uncertainty during migration. Note that you can use it only after obtaining and registering in you current session Google Api Key. For details on the API key check [here](http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error).
#'
#' @param Result FLightR result object obtained from \code{\link{run.particle.filter}}
#' @param dates Use NULL if all twilights will be used for plotting, one integer if specific twilight should be plotted (line number in Result$Results$Quantiles). Use data.frame with first column - start of the period and second - end of the period and each line represents a new period to plot specific periods, e.g. wintering or migration.
#' @param map.options options passed to \code{\link[ggmap]{get_map}}, note that \code{zoom} option is defined separately
#' @param percentiles Probability breaks for utilisation distribution
#' @param zoom Zoom for map. If 'auto' FLightR will try to find optimal zoom level by downloading different size maps and checking whether all the points fit the map.
#' @param geom_polygon.options options passed to \code{\link[ggplot2]{geom_polygon}}
#' @param save.options options passed to \code{\link[ggplot2]{ggsave}}. Filename should be defined here.
#' @param color.palette colors for probability contours. Either NULL or \code{\link[grDevices]{colorRampPalette}} object
#' @param use.palette should the same colors be used for polygon boundaries as for polygon filling?
#' @param background if provided will be used as a background. Must be created by \code{link[ggmap]{get_map}}
#' @param plot should function produce a plot?
#' @param save should function save results with \code{\link[ggplot2]{ggsave}}?
#' @return list with two parts
#' \item{res_buffers}{spatial buffers for defined probability values}
#' \item{p}{\code{\link[ggplot2]{ggplot}} object}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-06-25', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=1)
#' # here we will run only 1e4 partilces for a very short track.
#' # One should use 1e6 particles for the full run
#' Result<-run.particle.filter(all.in, threads=1,
#' nParticles=1e3, known.last=TRUE,
#' precision.sd=25, check.outliers=FALSE)
#'
#'\dontrun{
#' plot_util_distr(Result, zoom=6, save=FALSE)
#'}
#'
#' @author Eldar Rakhimberdiev
#' @export plot_util_distr
plot_util_distr<-function(Result, dates=NULL, map.options=NULL, percentiles=c(0.4, 0.6, 0.8), zoom="auto", geom_polygon.options=NULL, save.options=NULL, color.palette=NULL, use.palette=TRUE, background=NULL, plot=TRUE, save=TRUE) {
if (utils::packageVersion('ggmap')[1]<2.7) {stop('plot_util_distr function works only with ggmap >= 2.7.x')}
if (!ggmap::has_google_key()) stop('From August 2018 Google allows to use Google maps only for users with the API key, please get one and proceed as described here: http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error/')
Opt<-options('ggmap')
if(is.null(Opt$ggmap$google$account_type)) Opt$ggmap$google$account_type<-'standard'
if(is.null(Opt$ggmap$google$day_limit)) Opt$ggmap$google$day_limit<-2500
if(is.null(Opt$ggmap$google$second_limit)) Opt$ggmap$google$day_limit<-50
if(is.null(Opt$ggmap$google$client)) Opt$ggmap$google$client<-NA
if(is.null(Opt$ggmap$google$signature)) Opt$ggmap$google$signature<-NA
options('ggmap'= Opt$ggmap)
# for r cmd check
long<-NA
lat<-NA
group<-NA
## end ##########
if (is.null(color.palette)) color.palette <- grDevices::colorRampPalette(rev(c("#edf8fb",
"#b3cde3",
"#8c96c6",
"#88419d")))
# percentiles select area where bird spend at least xx percent of it's time
percentiles<-sort(percentiles)
r=get_gunion_r(Result)
#--------------
# Now I want to figure out minimum percentile plossible
res_buffers<-c()
nPoints<-c()
for (percentile in percentiles) {
res_cur<- get_time_spent_buffer(Result, dates, percentile, r)
res_buffers<-c(res_buffers, res_cur$Buffer)
nPoints<-c(nPoints, res_cur$nPoints)
}
# Now I want to take only the upper levels from identical nPoints
#Index_not_dupl<-(1:length(percentiles))[!duplicated(nPoints)]
#-----------------#
# and now we want to plot it
#if (is.null(background)) {
# combine location with map options
if (is.null(map.options)) map.options<-list()
#if (is.null(map.options$location)) map.options$location<-location
if (is.null(map.options$zoom) & is.numeric(zoom)) map.options$zoom=zoom
if (is.null(map.options$col)) map.options$col="bw"
#location<-res_buffers[[length(res_buffers)]]@bbox
location<-sf::st_bbox(res_buffers[[length(res_buffers)]])
#if (is.null(map.options$location)) map.options$location<-rowMeans(res_buffers[[length(res_buffers)]]@bbox)
if (is.null(map.options$location)) map.options$location<-sf::st_coordinates(sf::st_centroid(sf::st_as_sfc(sf::st_bbox(res_buffers[[length(res_buffers)]]))))
if (zoom=="auto") {
for (zoom_cur in (2:10)) {
map.options$zoom=zoom_cur
background <-do.call(ggmap::get_map, map.options)
bb<-attr(background, 'bb')
bb[2]<-ifelse(bb[2]< (-180), bb[2]+360, bb[2])
bb[4]<-ifelse(bb[4]< (-180), bb[4]+360, bb[4])
bb[2]<-ifelse(bb[2]> (180), bb[2]-360, bb[2])
bb[4]<-ifelse(bb[4]> (180), bb[4]-360, bb[4])
if (bb[4]<bb[2]) {
lonisinbb<-(location[1] >= bb[2] || location[1] <= bb[4]) & (location[3] >= bb[2] || location[3] <= bb[4])
} else {
lonisinbb<-(location[1] >= bb[2] && location[1] <= bb[4]) & (location[3] >= bb[2] && location[3] <= bb[4])
}
isinbb<-lonisinbb & location[2]>=bb[1] & location[4]<=bb[3]
if (!isinbb) break
}
map.options$zoom=zoom_cur-1
}
background <-do.call(ggmap::get_map, map.options)
#}
geom_polygon.options.external=geom_polygon.options
p<-ggmap::ggmap(background, maprange=TRUE)
#p<-ggmap::ggmap(background, extent = "normal", maprange=FALSE)
#p<-ggmap::ggmap(background, extent = "normal")
Colors<-color.palette(length(res_buffers))
for (i in length(res_buffers):1) {
buff_cur<-res_buffers[[i]]
if (is.null(geom_polygon.options)) geom_polygon.options=list()
if (is.null(geom_polygon.options$alpha)) geom_polygon.options$alpha=0.5
geom_path.options<-list()
geom_path.options$color<- geom_polygon.options$color
if (is.null(geom_path.options$color)) geom_path.options$color =ifelse(use.palette, Colors[i], '#756bb1')
if (is.null(geom_polygon.options$fill)) geom_polygon.options$fill =geom_path.options$color
geom_path.options$size<- geom_polygon.options$size
if (is.null(geom_path.options$size)) geom_path.options$size=2
geom_path.options$size<- geom_polygon.options$size
if (is.null(geom_path.options$size)) geom_path.options$size=2
geom_path.options$alpha=0.8
geom_path.options$linetype<- geom_polygon.options$linetype
geom_polygon.options$linetype=0
geom_polygon.options$mapping=ggplot2::aes(x=long, y=lat, group=group)
geom_path.options$mapping=geom_polygon.options$mapping
geom_polygon.options$data=buff_cur
geom_path.options$data=geom_polygon.options$data
p<-p+do.call(ggplot2::geom_polygon, geom_polygon.options)
p<-p+do.call(ggplot2::geom_path, geom_path.options)
geom_polygon.options=geom_polygon.options.external
}
if (plot) print(p)
if (save) {
if (is.null(save.options)) save.options=list()
if (is.null(save.options$filename)) save.options$filename<-"time_coverage.png"
save.options$plot <-p
save.options$dpi <-600
tmp<-do.call(ggplot2::ggsave, save.options)
}
b<-do.call(rbind, res_buffers)
SPDF = cbind(b, data.frame(percentile = percentiles, row.names=names(b)))
return(list(res_buffers=SPDF, p=p, bg=background))
}
seasonal_donut<-function() {
Seasonal_palette<-grDevices::colorRampPalette(grDevices::hsv(1-((1:365)+(365/4))%%365/365, s=0.8, v=0.8), space="Lab")
pie.data<-data.frame(group=as.factor(c(1:12)), value=rep(1,12))
pie.data$fraction = pie.data$value / sum(pie.data$value)
pie.data$ymax = cumsum(pie.data$fraction)
pie.data$ymin = c(0, utils::head(pie.data$ymax, n = -1))
donut<-ggplot2::ggplot(data = pie.data, ggplot2::aes(fill = pie.data$group, ymax = pie.data$ymax, ymin = pie.data$ymin, xmax = 1, xmin = 2)) +
ggplot2::geom_rect(colour = "grey30", show.legend = FALSE) +
ggplot2::coord_polar(theta = "y", start=pi) +
ggplot2::xlim(c(0, 2)) +
ggplot2::theme_bw() +
ggplot2::theme(plot.background = ggplot2::element_rect(fill = "transparent" ,colour = NA))+
ggplot2::theme(panel.background = ggplot2::element_blank())+
ggplot2::theme(panel.grid=ggplot2::element_blank()) +
ggplot2::theme(axis.text=ggplot2::element_blank()) +
ggplot2::theme(axis.ticks=ggplot2::element_blank()) +
ggplot2::theme(panel.border=ggplot2::element_blank()) +
ggplot2::scale_fill_manual(values=c(Seasonal_palette(13)))+
ggplot2::geom_text(ggplot2::aes(x = 1.5, y = (( pie.data$ymin+ pie.data$ymax)/2), label = pie.data$group), colour=grDevices::grey(0.99), size=5) +
ggplot2::xlab("") +
ggplot2::ylab("")
return(donut)
}
#' plot likelihood surface over map
#'
#' plots specific likelihood surface over map
#' @details function plots likelihoods before particle filter run, so these are pure results of calibrations without any movement model
#' @param object either output from \code{\link{make.prerun.object}} or \code{\link{run.particle.filter}}
#' @param date either NULL or a date (possibly with time) closest to the twilight you wan to be plotted
#' @param twilight.index number of likelihood surface to be plotted
#' @return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#' plot_likelihood(all.in, twilight.index=10)
#'
#' @author Eldar Rakhimberdiev
#' @export plot_likelihood
plot_likelihood<-function(object, date=NULL, twilight.index=NULL) {
my.golden.colors <- grDevices::colorRampPalette(c("white","#FF7100"))
if (!is.null(date)) {
date<-as.POSIXct(date, tz='GMT')
twilight.index<-which.min(abs(object$Indices$Matrix.Index.Table$time-date))
}
fields::image.plot(fields::as.image(object$Spatial$Phys.Mat[,twilight.index], x=object$Spatial$Grid[,1:2],nrow=60, ncol=60),
col=my.golden.colors(64), main=paste("twilight number",twilight.index, format(object$Indices$Matrix.Index.Table$time[twilight.index], tz='UTC')))
maps::map('world2', add=TRUE)
return(NULL)
}
get_points_distribution<-function(Result, twilights.index) {
Points_rle<-Result$Results$Points.rle[twilights.index]
All.Points<-rep(0, nrow(Result$Spatial$Grid))
for (twilight in 1:length(twilights.index)) {
All.Points[Points_rle[[twilight]]$values]<-All.Points[Points_rle[[twilight]]$values] + Points_rle[[twilight]]$lengths
}
return(All.Points)
}
get_utilisation_points<-function(Result, twilights.index, percentile) {
All.Points<-get_points_distribution(Result, twilights.index)
Order<-order(All.Points, decreasing=TRUE)
Order_rle<-rle(All.Points[Order]) #
nParticles<- sum(Result$Results$Points.rle[[1]]$lengths)
Cum_probs<-cumsum(Order_rle$values*Order_rle$length/nParticles/length(twilights.index))
tmp<-which(Cum_probs<=percentile)
if (length(tmp)==0) {
Selected_rle<-1 # if there is no such bin - return the most likely one
} else {
Selected_rle<-max(which(Cum_probs<=percentile))
}
Points<-Order[1:sum(Order_rle$lengths[1:Selected_rle])]
Points_selected<-(1:length(All.Points))[Points]
return(Points_selected)
}
eldarrak/FLightR documentation built on Nov. 10, 2023, 7:53 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions get_utilisation_points get_points_distribution plot_likelihood seasonal_donut plot_util_distr get_time_spent_buffer get_gunion_r get_buffer plot_lon_lat map.FLightR.ggmap
Documented in map.FLightR.ggmap plot_likelihood plot_lon_lat plot_util_distr
#' plots result over map
#'
#' plots track over map with probability cloud. Can plot only part of the track if dates are specified. Note that you can use it only after obtaining and registering in you current session Google Api Key. For details on the API key check [here](http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error).
#' @param Result FLightR result object obtained from \code{\link{run.particle.filter}}
#' @param dates either NULL if all twilights should be included or data.frame with first column - start of the period and second end of the period. Each line represents a new period
#' @param plot.cloud Should probability cloud be plotted? If TRUE cloud is estimated by \code{\link[ggplot2]{stat_density2d}}
#' @param map.options options passed to \code{\link[ggmap]{get_map}}, note that \code{zoom} option is defined separately
#' @param plot.options plotting options. Not defined yet!
#' @param save.options options passed to \code{\link[ggplot2]{ggsave}}. Filename should be defined here.
#' @param zoom Zoom for map. If 'auto' FLightR will try to find optimal zoom level by downloading different size maps and checking whether all the points fit the map.
#' @param return.ggobj Should ggobj be returned for subsequent checks and/or replotting
#' @param seasonal.colors if true points of the track will have seasonal colors
#' @param seasonal.donut.location if NULL - no color wheel placed, otherwise select one of 'bottomleft', 'bottomright', 'topleft'
#' @param seasonal.donut.proportion how much of X axis should color wheel occupy.
#' return either NULL or ggplot2 class object
#' @param save should function save results with \code{\link[ggplot2]{ggsave}}?
#' @return if 'return.ggobj=TRUE' return ggplot object otherwise returns 'NULL'.
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-06-25', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#' # here we will run only 1e4 partilces for a very short track.
#' # One should use 1e6 particles for the full run
#' Result<-run.particle.filter(all.in, threads=1,
#' nParticles=1e3, known.last=TRUE,
#' precision.sd=25, check.outliers=FALSE)
#'
#'\dontrun{
#' map.FLightR.ggmap(Result, seasonal.donut.location=NULL, zoom=6, save=FALSE)
#'}
#' # for this short track without variance seasonal donut does not work,
#' # but for normall track it will.
#' @author Eldar Rakhimberdiev
#' @export map.FLightR.ggmap
map.FLightR.ggmap<-function(Result, dates=NULL, plot.cloud=TRUE, map.options=NULL, plot.options=NULL, save.options=NULL, zoom="auto", return.ggobj=FALSE, seasonal.colors=TRUE, seasonal.donut.location='topleft', seasonal.donut.proportion=0.5, save=TRUE) {
if (!is.null(plot.options)) warning("plot options are not in use yet. Let me know what you would like to have here.")
if (utils::packageVersion('ggmap')[1]<2.7) { stop('map.FLightR.ggmap function works only with ggmap >= 2.7.x') }
if (!ggmap::has_google_key()) stop('From August 2018 Google allows to use Google maps only for users with the API key, please get one and proceed as described here: http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error/')
# dates should be a data.frame with first point - starting dates and last column end dates for periods
# ggsave.options is a list that will be will be directly passed to ggsave
Opt<-options('ggmap')
if(is.null(Opt$ggmap$google$account_type)) Opt$ggmap$google$account_type<-'standard'
if(is.null(Opt$ggmap$google$day_limit)) Opt$ggmap$google$day_limit<-2500
if(is.null(Opt$ggmap$google$second_limit)) Opt$ggmap$google$day_limit<-50
if(is.null(Opt$ggmap$google$client)) Opt$ggmap$google$client<-NA
if(is.null(Opt$ggmap$google$signature)) Opt$ggmap$google$signature<-NA
options('ggmap'= Opt$ggmap)
# select twilights to plot
if (is.null(dates)) {
twilights.index<-1:length(Result$Results$Points.rle)
} else {
twilights.index<-c()
for (segment in 1:nrow(dates)) {
twilights.index<-c(twilights.index, which(Result$Results$Quantiles$time>=dates[segment,1] & Result$Results$Quantiles$time<=dates[segment,2]))
}
if (length(twilights.index)==0) stop("dates do not overlap with the track time span!")
}
if (plot.cloud) {
Points_rle<-Result$Results$Points.rle[twilights.index]
All.Points<-rep(0, nrow(Result$Spatial$Grid))
for (twilight in 1:length(twilights.index)) {
All.Points[Points_rle[[twilight]]$values]<-All.Points[Points_rle[[twilight]]$values] + Points_rle[[twilight]]$lengths
}
Points<-Result$Spatial$Grid[All.Points>0,,drop=FALSE][sample.int(length(All.Points[All.Points>0]), size = 10000, replace = TRUE, prob = All.Points[All.Points>0]), 1:2]
}
# background map
# check whether Grid was over dateline:
overdateline<-ifelse(attr(Result$Spatial$Grid, 'left')> attr(Result$Spatial$Grid, 'right'), TRUE, FALSE)
if (overdateline) {
location<-cbind(
min(Result$Results$Quantiles$Medianlon[twilights.index][Result$Results$Quantiles$Medianlon[twilights.index]>0]),
min(Result$Results$Quantiles$Medianlat[twilights.index]),
max(Result$Results$Quantiles$Medianlon[twilights.index][Result$Results$Quantiles$Medianlon[twilights.index]<0]),
max(Result$Results$Quantiles$Medianlat[twilights.index]))
if (!is.finite(location[1])) location[1] <-180
if (!is.finite(location[3])) location[3] <- -180
} else {
location<-cbind(
min(Result$Results$Quantiles$Medianlon[twilights.index]),
min(Result$Results$Quantiles$Medianlat[twilights.index]),
max(Result$Results$Quantiles$Medianlon[twilights.index]),
max(Result$Results$Quantiles$Medianlat[twilights.index]))
}
# there could be problems here if the dateline is passed. I will add check for it later.
center_lon<-mean(c(location[1], location[3]))
if (overdateline) center_lon<-mean(c(location[1], location[3]+360))-360
if (center_lon < (-180)) center_lon<-center_lon+360
center_lat<-mean(c(location[2], location[4]))
if (center_lon>180 ) center_lon<-center_lon-360
if ((location[1])>180 ) location[1]<-location[1]-360
if ((location[3])>180 ) location[3]<-location[3]-360
# combine location with map options
if (is.null(map.options)) map.options<-list()
if (is.null(map.options$location)) map.options$location<-c(center_lon, center_lat)
if (is.null(map.options$zoom) & is.numeric(zoom)) map.options$zoom=zoom
if (is.null(map.options$col)) map.options$col="bw"
if (zoom=="auto") {
for (zoom_cur in (2:10)) {
map.options$zoom=zoom_cur
background <-do.call(ggmap::get_map, map.options)
bb<-attr(background, 'bb')
bb[2]<-ifelse(bb[2]< (-180), bb[2]+360, bb[2])
bb[4]<-ifelse(bb[4]< (-180), bb[4]+360, bb[4])
bb[2]<-ifelse(bb[2]> (180), bb[2]-360, bb[2])
bb[4]<-ifelse(bb[4]> (180), bb[4]-360, bb[4])
if (bb[4]<bb[2]) {
lonisinbb<-(location[1] >= bb[2] || location[1] <= bb[4]) & (location[3] >= bb[2] || location[3] <= bb[4])
} else {
lonisinbb<-(location[1] >= bb[2] && location[1] <= bb[4]) & (location[3] >= bb[2] && location[3] <= bb[4])
}
isinbb<-lonisinbb & location[2]>=bb[1] & location[4]<=bb[3]
if (!isinbb) break
}
map.options$zoom=zoom_cur-1
}
background <-do.call(ggmap::get_map, map.options)
p<-ggmap::ggmap(background)#, maprange=TRUE)
if (plot.cloud) {
#if (overdateline) Points[,1]<-ifelse(Points[,1]>0, Points[,1]-360, Points[,1])
# for r cmd check:
..level..<-NA
lon<-NA
lat<-NA
#################
p<- p+ ggplot2::stat_density2d(data=data.frame(Points), ggplot2::aes(fill = ..level.., alpha = ..level.., x=lon, y=lat), size = 0.01, geom = 'polygon', n=400)
p<-p+ ggplot2::scale_fill_gradient(low = "green", high = "red")
p<-p+ ggplot2::scale_alpha(range = c(0.00, 0.25), guide = FALSE)
}
#p<-p+ggplot2::coord_map(projection="mercator",
# xlim=c(attr(background, "bb")$ll.lon, attr(background, "bb")$ur.lon),
# ylim=c(attr(background, "bb")$ll.lat, attr(background, "bb")$ur.lat))
p<-p+ ggplot2::theme(legend.position = "none", axis.title = ggplot2::element_blank(), text = ggplot2::element_text(size = 12))
# here I plot track for selected dates
Points2plot<-data.frame(lat=Result$Results$Quantiles$Medianlat,lon=Result$Results$Quantiles$Medianlon)
if (overdateline) Points2plot[,2][Points2plot[,2]>attr(background, "bb")$ur.lon] <- Points2plot[,2][Points2plot[,2]>attr(background, "bb")$ur.lon]-360
if (overdateline) Points2plot[,2][Points2plot[,2]<attr(background, "bb")$ll.lon] <- Points2plot[,2][Points2plot[,2]<attr(background, "bb")$ll.lon] + 360
if (is.null(dates)) {
cur_twilights<-1:length(Result$Results$Quantiles$time)
Total_segments<-1
} else {
cur_twilights<-which(Result$Results$Quantiles$time>=dates[segment,1] & Result$Results$Quantiles$time<=dates[segment,2])
Total_segments<-nrow(dates)
}
for (segment in 1:Total_segments) {
p<-p+ ggplot2::geom_path(data=Points2plot[cur_twilights,],ggplot2::aes(x=lon,y=lat), colour=grDevices::grey(0.3))
if (seasonal.colors) {
Seasonal_palette<-grDevices::colorRampPalette(grDevices::hsv(1-((1:365)+(365/4))%%365/365, s=0.8, v=0.8), space="Lab")
Days<-as.integer(format(Result$Results$Quantiles$time[cur_twilights], "%j"))
p<-p+ggplot2::geom_point(data=Points2plot[cur_twilights,], shape=21, colour=grDevices::grey(0.3), bg=Seasonal_palette(365)[Days])
} else {
p<-p+ggplot2::geom_point(data=Points2plot[cur_twilights,], shape="+", colour=grDevices::grey(0.3))
}
}
if (!is.null(seasonal.donut.location)) {
Xrange<-c(attr(background, "bb")$ll.lon, attr(background, "bb")$ur.lon)
Yrange<-c(attr(background, "bb")$ll.lat, attr(background, "bb")$ur.lat)
d<-seasonal_donut()
g = ggplot2::ggplotGrob(d)
if (seasonal.donut.location=='bottomleft') {
p<-p + ggmap::inset(grob = g, xmin = Xrange[1], xmax = Xrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]), ymin = Yrange[1], ymax = Yrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]))
}
if (seasonal.donut.location=='topleft') {
p<-p + ggmap::inset(grob = g, xmin = Xrange[1], xmax = Xrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]), ymin = Yrange[2]-seasonal.donut.proportion*(Yrange[2]-Yrange[1]), ymax = Yrange[2])
}
if (seasonal.donut.location=='bottomright') {
p<-p + ggmap::inset(grob = g, xmin = Xrange[2]-seasonal.donut.proportion*(Yrange[2]-Yrange[1]), xmax = Xrange[2], ymin = Yrange[1], ymax = Yrange[1]+seasonal.donut.proportion*(Yrange[2]-Yrange[1]))
}
}
graphics::plot(p)
if (save) {
if (is.null(save.options)) save.options<-list()
if (is.null(save.options$filename)) save.options$filename<-"FLightR.map.pdf"
do.call(ggplot2::ggsave, save.options)
}
if (return.ggobj) {return(p)} else {return(NULL)}
}
#' plots result by longitude and latitude
#'
#' This function plots result by latitude and longitude in either vertical or horizontal layout.
#' @param Result FLightR result object obtained from \code{\link{run.particle.filter}}
#' @param scheme either 'vertical' or 'horizontal' layouts
#' return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#' # here we will run only 1e4 partilces for a very short track.
#' # One should use 1e6 particles for the full run
#' Result<-run.particle.filter(all.in, threads=1,
#' nParticles=1e3, known.last=TRUE,
#' precision.sd=25, check.outliers=FALSE)
#'
#' plot_lon_lat(Result)
#'
#' @author Eldar Rakhimberdiev
#' @export plot_lon_lat
plot_lon_lat<-function(Result, scheme=c("vertical", "horizontal")) {
oldpar <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(oldpar))
Quantiles<-Result$Results$Quantiles
# check whether Grid was over dateline:
overdateline<-ifelse(attr(Result$Spatial$Grid, 'left')> attr(Result$Spatial$Grid, 'right'), TRUE, FALSE)
if (overdateline) {
Quantiles$Medianlon[Quantiles$Medianlon<0]<-Quantiles$Medianlon[Quantiles$Medianlon<0]+360
Quantiles$LCI.lon[Quantiles$LCI.lon<0]<-Quantiles$LCI.lon[Quantiles$LCI.lon<0]+360
Quantiles$UCI.lon[Quantiles$UCI.lon<0]<-Quantiles$UCI.lon[Quantiles$UCI.lon<0]+360
Quantiles$TrdQu.lon[Quantiles$TrdQu.lon<0]<-Quantiles$TrdQu.lon[Quantiles$TrdQu.lon<0]+360
Quantiles$FstQu.lon[Quantiles$FstQu.lon<0]<-Quantiles$FstQu.lon[Quantiles$FstQu.lon<0]+360
}
if (scheme[1]=="horizontal") {
graphics::par(mfrow=c(1,2))
} else {
graphics::par(mfrow=c(2,1))
}
graphics::par(mar=c(2,4,3,1),cex=1)
suppressWarnings(Sys.setlocale("LC_ALL", "English"))
#------here I have create the equinox dates..
Years<-unique(format(Quantiles$time, format="%Y"))
eq<-c(as.POSIXct(paste(Years, "-09-22 00:00:00 GMT", sep="")), as.POSIXct(paste(Years, "-03-20 12:00:00 GMT", sep="")))
eq<-eq[eq>min(Quantiles$time) & eq<max(Quantiles$time)]
#------- vert grid
Vert_grid<-seq(as.POSIXct("2000-01-01", tz='GMT'), as.POSIXct("2030-01-01", tz='GMT'), by="month")
Vert_grid<-Vert_grid[Vert_grid>=min(Quantiles$time) & Vert_grid<=max(Quantiles$time)]
#Longitude
graphics::plot(Quantiles$Medianlon~Quantiles$time, las=1,col=grDevices::grey(0.1),pch=16,
ylab="Longitude",xlab="",lwd=2, ylim=range(c( Quantiles$LCI.lon,
Quantiles$UCI.lon )), type="n", axes=FALSE)
graphics::axis(2, las=1)
graphics::axis.POSIXct(1, x=Quantiles$time, format="1-%b")
graphics::box()
# add vertical lines for the first day of every month
graphics::abline(v=Vert_grid, col=grDevices::grey(0.5), lty=2)
graphics::abline(h=seq(-180, 360, by=10), col=grDevices::grey(0.5), lty=2)
graphics::abline(v=eq, col=2, lwd=2, lty=1)
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)),
y=c(Quantiles$LCI.lon, rev(Quantiles$UCI.lon)),
col=grDevices::grey(0.9), border=grDevices::grey(0.5))
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)),
y=c(Quantiles$TrdQu.lon, rev(Quantiles$FstQu.lon)),
col=grDevices::grey(0.7), border=grDevices::grey(0.5))
graphics::lines(Quantiles$Medianlon~Quantiles$time, col=grDevices::grey(0.1),lwd=2)
#Latitude
graphics::par(mar=c(3,4,1,1))
graphics::plot(Quantiles$Medianlat~Quantiles$time, las=1,col=grDevices::grey(0.1),
pch=16,ylab="Latitude",xlab="",lwd=2,
ylim=range(c( Quantiles$UCI.lat, Quantiles$LCI.lat )), type="n", axes=FALSE)
graphics::axis(2, las=1)
graphics::axis.POSIXct(1, x=Quantiles$time, format="1-%b")
graphics::box()
graphics::abline(v=Vert_grid, col=grDevices::grey(0.5), lty=2)
graphics::abline(h=seq(-80, 80, by=10), col=grDevices::grey(0.5), lty=2)
graphics::abline(v=eq, col=2, lwd=2, lty=1)
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)), y=c(Quantiles$LCI.lat, rev(Quantiles$UCI.lat)),
col=grDevices::grey(0.9), border=grDevices::grey(0.5))
graphics::polygon(x=c(Quantiles$time, rev(Quantiles$time)), y=c(Quantiles$TrdQu.lat, rev(Quantiles$FstQu.lat)),
col=grDevices::grey(0.7), border=grDevices::grey(0.5))
graphics::lines(Quantiles$Medianlat~Quantiles$time, col=grDevices::grey(0.1),lwd=2)
return(NULL)
}
get_buffer<-function(coords, r){
# p = sp::SpatialPoints(matrix(coords, ncol=2), proj4string=sp::CRS("+proj=longlat +datum=WGS84"))
# aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
# p@coords[[2]], p@coords[[1]])
# projected <- sp::spTransform(p, sp::CRS(aeqd))
# buffered <- rgeos::gBuffer(projected, width=r, byid=TRUE)
# buffer_lonlat <- sp::spTransform(buffered, CRS=p@proj4string)
p <- sf::st_as_sf(as.data.frame(matrix(coords, ncol=2)), coords=c('V1','V2'),crs = 4326)
aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
sf::st_coordinates(p)[[2]], sf::st_coordinates(p)[[1]])
projected <- sf::st_transform(p,sf::st_crs(aeqd))
buffered <- sf::st_buffer(projected, dist = r)
buffer_lonlat <- sf::st_transform(buffered,sf::st_crs(p))
return(buffer_lonlat)
}
get_gunion_r<-function(Result) {
#Distances= sp::spDists(Result$Spatial$Grid[1:min(c(nrow(Result$Spatial$Grid), 1000)),1:2], longlat=TRUE)
Distances= as.numeric(sf::st_distance(sf::st_as_sf(as.data.frame(Result$Spatial$Grid[1:min(c(nrow(Result$Spatial$Grid), 1000)),1:2]),
coords=c('lon','lat'),
crs=4326))/1000)
# ok, distances go up to 51.2.. the next step is 62..
# so if I round them
Selected_dist<-unique(sort(round(Distances/10)*10))[2]
r<-Selected_dist*0.85*1000
return(r)
}
get_time_spent_buffer<-function(Result, dates=NULL, percentile=0.5, r=NULL) {
# r in meters..
# dates could be either NULL - for plotting all available dates,
# or numeric with length of one - to plot a specific twilight by number (row number in Quantiles)
# or data.frame with two columns - one for the first date in period and one for the last date in period.
if (is.null(dates)) {
twilights.index<-1:length(Result$Results$Points.rle)
} else {
if (is.numeric(dates) & length(dates)==1) {
twilights.index<-dates
} else {
twilights.index<-c()
for (segment in 1:nrow(dates)) {
twilights.index<-c(twilights.index, which(Result$Results$Quantiles$time>=dates[segment,1] & Result$Results$Quantiles$time<=dates[segment,2]))
}
if (length(twilights.index)==0) stop("dates do not overlap with the track time span!")
}
}
message('function will plot ', length(twilights.index), ' twilights\n')
Points_selected<-get_utilisation_points(Result, twilights.index, percentile)
if (is.null(r)) {
r=get_gunion_r(Result)
}
#Now I want to create a spatial buffers around points with this radius
#spoints = SpatialPoints(Result$Spatial$Grid[Points_selected,1:2], proj4string= CRS("+proj=longlat +datum=WGS84"))
Buffers<-apply(matrix(Result$Spatial$Grid[Points_selected,1:2], ncol=2), 1, get_buffer, r=r)
Buff_comb<-Buffers[[1]]
if (length(Buffers)>1) {
for (i in 2:length(Buffers)) {
#Buff_comb<-rgeos::gUnion(Buff_comb, Buffers[[i]])
Buff_comb<-sf::st_union(Buff_comb, Buffers[[i]])
}
}
#Buff_comb_simpl<-rgeos::gSimplify(Buff_comb, tol=0.01, topologyPreserve=TRUE)
Buff_comb_simpl<-sf::st_simplify(Buff_comb,preserveTopology = FALSE, dTolerance = 0.01)
return(list(Buffer=Buff_comb_simpl, nPoints=length(Points_selected))) # Points before
}
#' plots resulting track over map with uncertainty shown by space utilisation distribution
#'
#' May be use not only for the whole track but for a set of specific dates, e.g. to show spatial uncertainty during migration. Note that you can use it only after obtaining and registering in you current session Google Api Key. For details on the API key check [here](http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error).
#'
#' @param Result FLightR result object obtained from \code{\link{run.particle.filter}}
#' @param dates Use NULL if all twilights will be used for plotting, one integer if specific twilight should be plotted (line number in Result$Results$Quantiles). Use data.frame with first column - start of the period and second - end of the period and each line represents a new period to plot specific periods, e.g. wintering or migration.
#' @param map.options options passed to \code{\link[ggmap]{get_map}}, note that \code{zoom} option is defined separately
#' @param percentiles Probability breaks for utilisation distribution
#' @param zoom Zoom for map. If 'auto' FLightR will try to find optimal zoom level by downloading different size maps and checking whether all the points fit the map.
#' @param geom_polygon.options options passed to \code{\link[ggplot2]{geom_polygon}}
#' @param save.options options passed to \code{\link[ggplot2]{ggsave}}. Filename should be defined here.
#' @param color.palette colors for probability contours. Either NULL or \code{\link[grDevices]{colorRampPalette}} object
#' @param use.palette should the same colors be used for polygon boundaries as for polygon filling?
#' @param background if provided will be used as a background. Must be created by \code{link[ggmap]{get_map}}
#' @param plot should function produce a plot?
#' @param save should function save results with \code{\link[ggplot2]{ggsave}}?
#' @return list with two parts
#' \item{res_buffers}{spatial buffers for defined probability values}
#' \item{p}{\code{\link[ggplot2]{ggplot}} object}
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-06-25', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=1)
#' # here we will run only 1e4 partilces for a very short track.
#' # One should use 1e6 particles for the full run
#' Result<-run.particle.filter(all.in, threads=1,
#' nParticles=1e3, known.last=TRUE,
#' precision.sd=25, check.outliers=FALSE)
#'
#'\dontrun{
#' plot_util_distr(Result, zoom=6, save=FALSE)
#'}
#'
#' @author Eldar Rakhimberdiev
#' @export plot_util_distr
plot_util_distr<-function(Result, dates=NULL, map.options=NULL, percentiles=c(0.4, 0.6, 0.8), zoom="auto", geom_polygon.options=NULL, save.options=NULL, color.palette=NULL, use.palette=TRUE, background=NULL, plot=TRUE, save=TRUE) {
if (utils::packageVersion('ggmap')[1]<2.7) {stop('plot_util_distr function works only with ggmap >= 2.7.x')}
if (!ggmap::has_google_key()) stop('From August 2018 Google allows to use Google maps only for users with the API key, please get one and proceed as described here: http://ornithologyexchange.org/forums/topic/38315-mapflightrggmap-error/')
Opt<-options('ggmap')
if(is.null(Opt$ggmap$google$account_type)) Opt$ggmap$google$account_type<-'standard'
if(is.null(Opt$ggmap$google$day_limit)) Opt$ggmap$google$day_limit<-2500
if(is.null(Opt$ggmap$google$second_limit)) Opt$ggmap$google$day_limit<-50
if(is.null(Opt$ggmap$google$client)) Opt$ggmap$google$client<-NA
if(is.null(Opt$ggmap$google$signature)) Opt$ggmap$google$signature<-NA
options('ggmap'= Opt$ggmap)
# for r cmd check
long<-NA
lat<-NA
group<-NA
## end ##########
if (is.null(color.palette)) color.palette <- grDevices::colorRampPalette(rev(c("#edf8fb",
"#b3cde3",
"#8c96c6",
"#88419d")))
# percentiles select area where bird spend at least xx percent of it's time
percentiles<-sort(percentiles)
r=get_gunion_r(Result)
#--------------
# Now I want to figure out minimum percentile plossible
res_buffers<-c()
nPoints<-c()
for (percentile in percentiles) {
res_cur<- get_time_spent_buffer(Result, dates, percentile, r)
res_buffers<-c(res_buffers, res_cur$Buffer)
nPoints<-c(nPoints, res_cur$nPoints)
}
# Now I want to take only the upper levels from identical nPoints
#Index_not_dupl<-(1:length(percentiles))[!duplicated(nPoints)]
#-----------------#
# and now we want to plot it
#if (is.null(background)) {
# combine location with map options
if (is.null(map.options)) map.options<-list()
#if (is.null(map.options$location)) map.options$location<-location
if (is.null(map.options$zoom) & is.numeric(zoom)) map.options$zoom=zoom
if (is.null(map.options$col)) map.options$col="bw"
#location<-res_buffers[[length(res_buffers)]]@bbox
location<-sf::st_bbox(res_buffers[[length(res_buffers)]])
#if (is.null(map.options$location)) map.options$location<-rowMeans(res_buffers[[length(res_buffers)]]@bbox)
if (is.null(map.options$location)) map.options$location<-sf::st_coordinates(sf::st_centroid(sf::st_as_sfc(sf::st_bbox(res_buffers[[length(res_buffers)]]))))
if (zoom=="auto") {
for (zoom_cur in (2:10)) {
map.options$zoom=zoom_cur
background <-do.call(ggmap::get_map, map.options)
bb<-attr(background, 'bb')
bb[2]<-ifelse(bb[2]< (-180), bb[2]+360, bb[2])
bb[4]<-ifelse(bb[4]< (-180), bb[4]+360, bb[4])
bb[2]<-ifelse(bb[2]> (180), bb[2]-360, bb[2])
bb[4]<-ifelse(bb[4]> (180), bb[4]-360, bb[4])
if (bb[4]<bb[2]) {
lonisinbb<-(location[1] >= bb[2] || location[1] <= bb[4]) & (location[3] >= bb[2] || location[3] <= bb[4])
} else {
lonisinbb<-(location[1] >= bb[2] && location[1] <= bb[4]) & (location[3] >= bb[2] && location[3] <= bb[4])
}
isinbb<-lonisinbb & location[2]>=bb[1] & location[4]<=bb[3]
if (!isinbb) break
}
map.options$zoom=zoom_cur-1
}
background <-do.call(ggmap::get_map, map.options)
#}
geom_polygon.options.external=geom_polygon.options
p<-ggmap::ggmap(background, maprange=TRUE)
#p<-ggmap::ggmap(background, extent = "normal", maprange=FALSE)
#p<-ggmap::ggmap(background, extent = "normal")
Colors<-color.palette(length(res_buffers))
for (i in length(res_buffers):1) {
buff_cur<-res_buffers[[i]]
if (is.null(geom_polygon.options)) geom_polygon.options=list()
if (is.null(geom_polygon.options$alpha)) geom_polygon.options$alpha=0.5
geom_path.options<-list()
geom_path.options$color<- geom_polygon.options$color
if (is.null(geom_path.options$color)) geom_path.options$color =ifelse(use.palette, Colors[i], '#756bb1')
if (is.null(geom_polygon.options$fill)) geom_polygon.options$fill =geom_path.options$color
geom_path.options$size<- geom_polygon.options$size
if (is.null(geom_path.options$size)) geom_path.options$size=2
geom_path.options$size<- geom_polygon.options$size
if (is.null(geom_path.options$size)) geom_path.options$size=2
geom_path.options$alpha=0.8
geom_path.options$linetype<- geom_polygon.options$linetype
geom_polygon.options$linetype=0
geom_polygon.options$mapping=ggplot2::aes(x=long, y=lat, group=group)
geom_path.options$mapping=geom_polygon.options$mapping
geom_polygon.options$data=buff_cur
geom_path.options$data=geom_polygon.options$data
p<-p+do.call(ggplot2::geom_polygon, geom_polygon.options)
p<-p+do.call(ggplot2::geom_path, geom_path.options)
geom_polygon.options=geom_polygon.options.external
}
if (plot) print(p)
if (save) {
if (is.null(save.options)) save.options=list()
if (is.null(save.options$filename)) save.options$filename<-"time_coverage.png"
save.options$plot <-p
save.options$dpi <-600
tmp<-do.call(ggplot2::ggsave, save.options)
}
b<-do.call(rbind, res_buffers)
SPDF = cbind(b, data.frame(percentile = percentiles, row.names=names(b)))
return(list(res_buffers=SPDF, p=p, bg=background))
}
seasonal_donut<-function() {
Seasonal_palette<-grDevices::colorRampPalette(grDevices::hsv(1-((1:365)+(365/4))%%365/365, s=0.8, v=0.8), space="Lab")
pie.data<-data.frame(group=as.factor(c(1:12)), value=rep(1,12))
pie.data$fraction = pie.data$value / sum(pie.data$value)
pie.data$ymax = cumsum(pie.data$fraction)
pie.data$ymin = c(0, utils::head(pie.data$ymax, n = -1))
donut<-ggplot2::ggplot(data = pie.data, ggplot2::aes(fill = pie.data$group, ymax = pie.data$ymax, ymin = pie.data$ymin, xmax = 1, xmin = 2)) +
ggplot2::geom_rect(colour = "grey30", show.legend = FALSE) +
ggplot2::coord_polar(theta = "y", start=pi) +
ggplot2::xlim(c(0, 2)) +
ggplot2::theme_bw() +
ggplot2::theme(plot.background = ggplot2::element_rect(fill = "transparent" ,colour = NA))+
ggplot2::theme(panel.background = ggplot2::element_blank())+
ggplot2::theme(panel.grid=ggplot2::element_blank()) +
ggplot2::theme(axis.text=ggplot2::element_blank()) +
ggplot2::theme(axis.ticks=ggplot2::element_blank()) +
ggplot2::theme(panel.border=ggplot2::element_blank()) +
ggplot2::scale_fill_manual(values=c(Seasonal_palette(13)))+
ggplot2::geom_text(ggplot2::aes(x = 1.5, y = (( pie.data$ymin+ pie.data$ymax)/2), label = pie.data$group), colour=grDevices::grey(0.99), size=5) +
ggplot2::xlab("") +
ggplot2::ylab("")
return(donut)
}
#' plot likelihood surface over map
#'
#' plots specific likelihood surface over map
#' @details function plots likelihoods before particle filter run, so these are pure results of calibrations without any movement model
#' @param object either output from \code{\link{make.prerun.object}} or \code{\link{run.particle.filter}}
#' @param date either NULL or a date (possibly with time) closest to the twilight you wan to be plotted
#' @param twilight.index number of likelihood surface to be plotted
#' @return 'NULL'
#' @examples
#' File<-system.file("extdata", "Godwit_TAGS_format.csv", package = "FLightR")
#' # to run example fast we will cut the real data file by 2013 Aug 20
#' Proc.data<-get.tags.data(File, end.date=as.POSIXct('2013-07-02', tz='GMT'))
#' Calibration.periods<-data.frame(
#' calibration.start=as.POSIXct(c(NA, "2014-05-05"), tz='GMT'),
#' calibration.stop=as.POSIXct(c("2013-08-20", NA), tz='GMT'),
#' lon=5.43, lat=52.93)
#' #use c() also for the geographic coordinates, if you have more than one calibration location
#' # (e. g., lon=c(5.43, 6.00), lat=c(52.93,52.94))
#'
#' # NB Below likelihood.correction is set to FALSE for fast run!
#' # Leave it as default TRUE for real examples
#' Calibration<-make.calibration(Proc.data, Calibration.periods, likelihood.correction=FALSE)
#'
#' Grid<-make.grid(left=0, bottom=50, right=10, top=56,
#' distance.from.land.allowed.to.use=c(-Inf, Inf),
#' distance.from.land.allowed.to.stay=c(-Inf, Inf))
#'
#' all.in<-make.prerun.object(Proc.data, Grid, start=c(5.43, 52.93),
#' Calibration=Calibration, threads=2)
#' plot_likelihood(all.in, twilight.index=10)
#'
#' @author Eldar Rakhimberdiev
#' @export plot_likelihood
plot_likelihood<-function(object, date=NULL, twilight.index=NULL) {
my.golden.colors <- grDevices::colorRampPalette(c("white","#FF7100"))
if (!is.null(date)) {
date<-as.POSIXct(date, tz='GMT')
twilight.index<-which.min(abs(object$Indices$Matrix.Index.Table$time-date))
}
fields::image.plot(fields::as.image(object$Spatial$Phys.Mat[,twilight.index], x=object$Spatial$Grid[,1:2],nrow=60, ncol=60),
col=my.golden.colors(64), main=paste("twilight number",twilight.index, format(object$Indices$Matrix.Index.Table$time[twilight.index], tz='UTC')))
maps::map('world2', add=TRUE)
return(NULL)
}
get_points_distribution<-function(Result, twilights.index) {
Points_rle<-Result$Results$Points.rle[twilights.index]
All.Points<-rep(0, nrow(Result$Spatial$Grid))
for (twilight in 1:length(twilights.index)) {
All.Points[Points_rle[[twilight]]$values]<-All.Points[Points_rle[[twilight]]$values] + Points_rle[[twilight]]$lengths
}
return(All.Points)
}
get_utilisation_points<-function(Result, twilights.index, percentile) {
All.Points<-get_points_distribution(Result, twilights.index)
Order<-order(All.Points, decreasing=TRUE)
Order_rle<-rle(All.Points[Order]) #
nParticles<- sum(Result$Results$Points.rle[[1]]$lengths)
Cum_probs<-cumsum(Order_rle$values*Order_rle$length/nParticles/length(twilights.index))
tmp<-which(Cum_probs<=percentile)
if (length(tmp)==0) {
Selected_rle<-1 # if there is no such bin - return the most likely one
} else {
Selected_rle<-max(which(Cum_probs<=percentile))
}
Points<-Order[1:sum(Order_rle$lengths[1:Selected_rle])]
Points_selected<-(1:length(All.Points))[Points]
return(Points_selected)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.