Merge pull request #113 from pedro-andrade-inpe/master

new functionality and small updates

DESCRIPTION

@@ -28,9 +28,12 @@ Imports:
       Rcpp,
       units,
       sf,
+      sp,
+      rgdal,
       sfheaders,
       lwgeom,
       utils,
+      raster,
       pbapply
 RoxygenNote: 7.0.2
 VignetteBuilder: knitr

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(append_height)
 export(cpp_snap_points)
 export(filter_by_agency_id)
 export(filter_by_shape_id)
@@ -19,7 +20,7 @@ importFrom(Rcpp,compileAttributes)
 importFrom(data.table,"%between%")
 importFrom(data.table,":=")
 importFrom(data.table,fifelse)
-importFrom(lwgeom,st_make_valid)
+importFrom(lwgeom,st_geod_length)
 importFrom(magrittr,"%>%")
 importFrom(stats,na.omit)
 importFrom(utils,head)

R/append_height.R

@@ -0,0 +1,27 @@
+#' @title Add a column with height to GPS data
+#' @description Add a column named height to GPS data using a tif data as reference.
+#' @param gps A GPS data created from gtfs2gps().
+#' @param heightfile The pathname of a tif file with height data.
+#' @return The GPS data with a new column named height.
+#' @export
+#' @examples
+#' library(dplyr)
+#' 
+#' fortaleza <- system.file("extdata/fortaleza.zip", package = "gtfs2gps")
+#' srtmfile <- system.file("extdata/fortaleza-srtm.tif", package = "gtfs2gps")
+#' 
+#' gtfs <- read_gtfs(fortaleza) %>%
+#'   filter_week_days() %>%
+#'   filter_single_trip()
+#' 
+#' fortaleza_gps <- gtfs2gps(gtfs, progress = FALSE) %>% append_height(srtmfile)
+append_height <- function(gps, heightfile){
+  f_gps <- gps
+
+  sp::coordinates(f_gps) <- ~shape_pt_lon + shape_pt_lat
+
+  myraster <- raster::raster(heightfile)
+  gps$height <- raster::extract(myraster, f_gps)
+
+  return(gps)
+}

R/gps_as_sf.R

@@ -2,8 +2,8 @@
 #' @title Convert GPS-like data.table to a Simple Feature object
 #' @description Convert a GPS data stored in a data.table into a Simple Feature.
 #' @param gps A data.table with timestamp data.
-#' @param crs A Coordinate Reference System. The default value is 4326 (latlong WGS84)
-#' @return A simple feature (sf) object.
+#' @param crs A Coordinate Reference System. The default value is 4326 (latlong WGS84).
+#' @return A simple feature (sf) object with point data.
 #' @export
 #' @examples
 #' library(dplyr)
@@ -16,10 +16,12 @@
 #' poa_gps_sf <- gps_as_sf(poa_gps)
 gps_as_sf <- function(gps, crs = 4326){
   # convert to sf
-  temp_gps <- sfheaders::sf_multipoint(gps, x = "shape_pt_lon", y = "shape_pt_lat",
-                                       multipoint_id = "shape_id", keep = TRUE)
+  temp_gps <- sfheaders::sf_point(gps, x = "shape_pt_lon", y = "shape_pt_lat", keep = TRUE)
 
-  temp_gps<- temp_gps[, -duplicated(names(temp_gps))]
+  dup <- duplicated(names(temp_gps))
+
+  if(any(dup))
+    temp_gps<- temp_gps[, -dup]
 
   # add projection
   sf::st_crs(temp_gps) <- crs

R/zzz.R

@@ -16,7 +16,7 @@ utils::globalVariables(c(".", "%>%", ":="))
 #' @importFrom utils head tail object.size
 #' @importFrom stats na.omit
 #' @importFrom Rcpp compileAttributes
-#' @importFrom lwgeom st_make_valid
+#' @importFrom lwgeom st_geod_length
 #' @useDynLib gtfs2gps, .registration = TRUE
 NULL
 

demo/fortaleza.R

@@ -6,19 +6,38 @@ fortaleza <- read_gtfs(system.file("extdata/fortaleza.zip", package="gtfs2gps"))
 names(fortaleza)
 
 ids <- fortaleza$shapes$shape_id %>% unique()
+ids
 
-fortaleza <- filter_by_shape_id(fortaleza, ids[1:2])
+fortaleza <- filter_by_shape_id(fortaleza, ids[1:4])
 
 for_small_sf <- gtfs_shapes_as_sf(fortaleza)
 
-for_gps <- gtfs2gps(fortaleza, parallel = TRUE)
+srtmfile <- system.file("extdata/fortaleza-srtm.tif", package = "gtfs2gps")
+
+for_gps <- gtfs2gps(fortaleza, parallel = TRUE) %>% append_height(srtmfile)
 
 for_gps_sf <- gps_as_sf(for_gps)
-for_gps_small <- for_gps[1:60, ]
+
+#write_sf(for_gps_sf, "for-gps2.shp")
+
+for_gps_small <- for_gps
 
 for_gps_small_sf <- gps_as_sf(for_gps_small)
 
 plot(sf::st_geometry(for_gps_small_sf), pch = 20)
 plot(sf::st_geometry(for_small_sf), col = "blue", add = TRUE)
 plot(sf::st_geometry(for_gps_sf), add = TRUE, pch = 20)
 box()
+
+# plotting GPS points with height
+
+srtmraster <- raster::raster(srtmfile)
+
+XYData<- data.frame(
+  X = for_gps$shape_pt_lon,
+  Y = for_gps$shape_pt_lat)
+
+plot(srtmraster)
+points(XYData, type = "p")
+
+plot(for_gps_sf["height"], pch = 20)

prep_data/fortaleza-srtm/generate-fortaleza-srtm-simplified.R

@@ -0,0 +1,18 @@
+# SRTM data downloaded from
+# http://www.dsr.inpe.br/topodata/acesso.php
+
+require(raster)
+require(gtfs2gps)
+
+inputDir <- "prep_data/fortaleza-srtm/"
+
+myraster <- paste0(inputDir, "03S39_ZN.tif") %>% raster::raster()
+
+fortaleza <- read_gtfs(system.file("extdata/fortaleza.zip", package="gtfs2gps"))
+
+fortaleza_sf <- gtfs_shapes_as_sf(fortaleza) %>% sf::st_buffer(0.001) %>% sf::as_Spatial()
+
+mm.sub <- raster::crop(myraster, fortaleza_sf)
+
+raster::writeRaster(mm.sub, "inst/extdata/fortaleza-srtm.tif", overwrite = TRUE)
+

prep_data/fortaleza.R

@@ -12,7 +12,7 @@ gtfs_shapes_as_sf(gtfs) %>%
 
 unique(gtfs$shapes$shape_id)
 
-shape_ids <- paste0("shape", 800:900, "-I")
+shape_ids <- paste0("shape", c(800:822, 825:854, 856:900), "-I")
 
 new_gtfs <- filter_by_shape_id(gtfs, shape_ids)
 

tests/testthat/test_append_height.R

@@ -0,0 +1,18 @@
+
+test_that("append_height", {
+  fortaleza <- system.file("extdata/fortaleza.zip", package="gtfs2gps")
+  srtmfile <- system.file("extdata/fortaleza-srtm.tif", package="gtfs2gps")
+
+  gtfs <- read_gtfs(fortaleza) %>%
+    filter_week_days() %>%
+    filter_single_trip()
+
+  fortaleza_gps <- gtfs2gps(gtfs, progress = FALSE) %>% append_height(srtmfile)
+
+  expect_equal(sum(fortaleza_gps$height), 72932, 0.05)
+
+  fort <- fortaleza_gps[which(is.na(fortaleza_gps$height)),]
+
+  #fott <- gps_as_sf(fort)
+  #sf::write_sf(fott, "fort-na.shp")  
+})