Merge pull request #27 from eldarrak/move2sf

Move2sf

.github/workflows/test-coverage.yaml

@@ -2,9 +2,9 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, master]
+    branches: [main, master, Move2sf]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, Move2sf]
 
 name: test-coverage
 

DESCRIPTION

@@ -2,12 +2,13 @@ Package: FLightR
 Type: Package
 Title: Reconstruct Animal Paths from Solar Geolocation Loggers Data
 Version: 0.5.4
-Date: 2023-09-xx
+Date: 2023-09-26
 Authors@R: c(
     person("Eldar", "Rakhimberdiev", email = "eldar.rakhimberdiev@uva.nl", role = c("aut","cre")),
     person("Anatoly", "Saveliev", role = "aut"),
     person("Julia", "Karagicheva", role = "aut"),
-    person("Simeon", "Lisovski", role = "ctb"))
+    person("Simeon", "Lisovski", role = "ctb"),
+	person("Johannes", "de Groeve", email = "j.degroeve@uva.nl", role = "ctb"))
 Description: Spatio-temporal locations of an animal are computed 
     from annotated data with a hidden Markov  model via particle
     filter algorithm. The package is relatively robust to varying
@@ -18,24 +19,24 @@ Description: Spatio-temporal locations of an animal are computed
 URL: https://CRAN.R-project.org/package=FLightR
 BugReports: https://github.com/eldarrak/FLightR/issues
 Depends:
-    R (>= 3.0.2)
+    R (>= 4.1.0)
 Imports:
     bit,
+	geosphere,
     ggsn,
     ggmap,
     ggplot2,
     CircStats,
     circular,
     fields,
-    maptools,
+	maps,
     methods,
     mgcv,
     nlme,
     parallel,
     RcppArmadillo,
-    rgdal,
-    rgeos,
-    sp,
+	sf,
+	suntools,
     truncnorm
 License: GPL-3
 ByteCompile: true

R/data_preparation.R

@@ -183,7 +183,7 @@ make.calibration<-function(Proc.data, Calibration.periods, model.ageing=FALSE, p
 find.stationary.location<-function(Proc.data, calibration.start,  calibration.stop, plot=TRUE, initial.coords=NULL, print.optimization=TRUE, reltol=1e-4) {
    if (is.null(initial.coords)) stop('current function vesrion requires some inital coordinates to start search, they should not be very close but within few thousand km!')
    ll_function<-function(initial.coords, Proc.data, calibration.start, calibration.stop, plot=TRUE, stage=1) {
-   sink("tmp")
+   sink("flightr_tmp")
         Calibration.period<-data.frame(
         calibration.start=as.POSIXct(calibration.start, tz='GMT'),
         calibration.stop=as.POSIXct(calibration.stop, tz='GMT'),
@@ -827,11 +827,11 @@ All.Days.extended<-seq(min(All.Days)-86400, max(All.Days)+86400, by="days")
    start_no_polar[2]<-min(abs(start_no_polar[2]), 65.73) * sign(start_no_polar[2])
 
 # these are all potential twilights that we would have for the start point...
-Potential.twilights<-sort(c(maptools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunrise", POSIXct.out=TRUE)[,2], maptools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunset", POSIXct.out=TRUE)[,2]))
+Potential.twilights<-sort(c(suntools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunrise", POSIXct.out=TRUE)[,2], suntools::sunriset(matrix(start_no_polar, nrow=1), All.Days.extended, direction="sunset", POSIXct.out=TRUE)[,2]))
 # now we need to add dask or dawn to it..
-Sunrise<-maptools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunrise", POSIXct.out=TRUE)[,2]
+Sunrise<-suntools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunrise", POSIXct.out=TRUE)[,2]
 Sunrises<-c(Sunrise-(3600*24), Sunrise, Sunrise+(3600*24))
-Sunset<-maptools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunset", POSIXct.out=TRUE)[,2]
+Sunset<-suntools::sunriset(matrix(start_no_polar, nrow=1), Potential.twilights[1], direction="sunset", POSIXct.out=TRUE)[,2]
 Sunsets<-c(Sunset-(3600*24), Sunset, Sunset+(3600*24))
 First.twilight<-ifelse(which.min(c(min(abs(difftime(Sunrises,Potential.twilights[1], units="mins"))), min(abs(difftime(Sunsets,Potential.twilights[1], units="mins")))))==1, "dawn", "dusk")
 
@@ -868,7 +868,10 @@ for (i in 1:length(processed.light$Final.dawn$Data$gmt)) {
 while (is.na(Index.tab$Curr.mat[1])) Index.tab<-Index.tab[-1,]
 while (is.na(Index.tab$Curr.mat[nrow(Index.tab)])) Index.tab<-Index.tab[-nrow(Index.tab),]
 Index.tab$Point<-NA
-First.Point<-which.min(sp::spDistsN1(Grid[,1:2], start,  longlat=TRUE))
+#First.Point<-which.min(sp::spDistsN1(Grid[,1:2], start,  longlat=TRUE))
+Grid_sf<-Grid[,1:2] |> sf::st_multipoint() |>  sf::st_sfc() |>   sf::st_cast('POINT')
+First.Point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
+
 Index.tab$Point[1]<-First.Point
 #
 # I decided that Curr.mat is not needed anymore
@@ -930,11 +933,15 @@ geologger.sampler.create.arrays<-function(Index.tab, Grid, start, stop=start) {
 
 	output$Spatial$Behav.mask<-as.integer(Grid[,3])
 
-	output$Spatial$start.point<-which.min(sp::spDistsN1(Grid[,1:2], start,  longlat=TRUE))
-
+	#output$Spatial$start.point<-which.min(sp::spDistsN1(Grid[,1:2], start,  longlat=TRUE))
+    Grid_sf<-Grid[,1:2] |> sf::st_multipoint() |>  sf::st_sfc() |>   sf::st_cast('POINT')
+    output$Spatial$start.point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
+
 	output$Spatial$start.location<-start
 	if (!is.na(stop[1]) ) {
-	   output$Spatial$stop.point<-which.min(sp::spDistsN1(Grid[,1:2], stop,  longlat=TRUE))
+	   #output$Spatial$stop.point<-which.min(sp::spDistsN1(Grid[,1:2], stop,  longlat=TRUE))
+	   output$Spatial$stop.point<- output$Spatial$start.point<-sf::st_nearest_feature(sf::st_point(start),Grid_sf, longlat = TRUE)
+
 	   output$Spatial$stop.location<-stop
 	} else {
 	   output$Spatial$stop.point<-NA

R/new_plotting_functions.R

@@ -349,17 +349,28 @@ plot_lon_lat<-function(Result, scheme=c("vertical", "horizontal")) {
 
 
 get_buffer<-function(coords, r){
-  p = sp::SpatialPoints(matrix(coords, ncol=2), proj4string=sp::CRS("+proj=longlat +datum=WGS84"))
+  # 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",
-                    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)
+                  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=  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]
@@ -400,11 +411,13 @@ get_time_spent_buffer<-function(Result, dates=NULL, percentile=0.5, r=NULL) {
    Buff_comb<-Buffers[[1]] 
    if (length(Buffers)>1) {
    for (i in 2:length(Buffers)) {
-       Buff_comb<-rgeos::gUnion(Buff_comb, Buffers[[i]])
+       #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)
-  return(list(Buffer=Buff_comb_simpl, nPoints=length(Points)))
+  #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
@@ -517,10 +530,12 @@ for (percentile in percentiles) {
 	#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
-
-    if (is.null(map.options$location)) map.options$location<-rowMeans(res_buffers[[length(res_buffers)]]@bbox)
-
+	#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)) {
@@ -618,29 +633,10 @@ for (percentile in percentiles) {
 	   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)))
 
-    my.rbind.SpatialPolygons = function(..., makeUniqueIDs = FALSE) {
-       dots = list(...)
-       names(dots) <- NULL
-       stopifnot(sp::identicalCRS(dots))
-       # checkIDSclash(dots)
-       pl = do.call(c, lapply(dots, function(x) methods::slot(x, "polygons")))
-       if (makeUniqueIDs)
-               pl = makeUniqueIDs(pl)
-       sp::SpatialPolygons(pl, proj4string = sp::CRS(sp::proj4string(dots[[1]])))
-	}
-    makeUniqueIDs <- function(lst) {
-	   ids = sapply(lst, function(i) methods::slot(i, "ID"))
-	   if (any(duplicated(ids))) {
-		  ids <- make.unique(as.character(unlist(ids)), sep = "")
-		  for (i in seq(along = ids))
-			lst[[i]]@ID = ids[i]
-	   }
-	   lst
-    }
-
-	b<-do.call(my.rbind.SpatialPolygons,  c(res_buffers, list(makeUniqueIDs=TRUE))) 
-    SPDF = sp::SpatialPolygonsDataFrame(b, data.frame(percentile = percentiles, row.names=names(b)))
 	return(list(res_buffers=SPDF, p=p, bg=background))
 }
 
@@ -714,9 +710,7 @@ plot_likelihood<-function(object, date=NULL, twilight.index=NULL) {
 
    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')))     
-   wrld_simpl<-NA				   
-   load(system.file("data", "wrld_simpl.rda", package = "maptools"))
-   sp::plot(wrld_simpl, add=TRUE)
+   maps::map('world2', add=TRUE)
    return(NULL)
 }