rdplyrgeospatialr-sfterra

Issue with setting CRS with make_track() in amt


I am having trouble setting the CRS with the crs argument within the make_track() function. It consistently just returns CRS as NA with no error. I am not sure if there is conflicting packages or if I am missing something.

library(sf)
library(terra)
library(amt) 
library(tibble)
library(dplyr)

# Load data from working directory, previously unzipped and downloaded from
# https://alaska.usgs.gov/data1/polarBear/polarBear_onLand_westHudsonBay_pagano/polarBear_GPSLocation_westHudsonBay_2019-2022.zip
polarBear_GPS <- read.csv("/Users/elleryvincent/Desktop/Spatial Ecology/polarBear_GPSLocation_westHudsonBay_2019-2022/polarBear_GPSLocation_westHudsonBay_2019-2022.csv")

# Check for NA values, 0 == no NA values in df
sum(is.na(polarBear_GPS))
# [1] 0

### Skip following if previous step returns 0

# ind <- complete.cases(polarBear_GPS[, c("Datetime", "Longitude", "Latitude")])
# 
# polarBear_GPS <- polarBear_GPS[ind == TRUE, ]
# 
# # Remove NA values from Datetime
# polarBear_GPS <- polarBear_GPS[!is.na(polarBear_GPS$Datetime), ]

###

# Make Datetime a date/time variable
polarBear_GPS$Datetime <- as.POSIXct(polarBear_GPS$Datetime, format="%m/%d/%Y %H:%M:%S", tz = "UTC")

# Check for duplicated timestamps per bear
polarBear_GPS |>
  summarise(ts_dupes = any(duplicated(Datetime)), .by = Bear)
 
# If any return true, run following step, else skip
polarBear_GPS <- polarBear_GPS |>
  distinct(Datetime, .keep_all = TRUE, .by = Bear)
st_crs(polarBear_GPS, 4326)

# Create the track with correct CRS (see metadata at download link above)
trk1 <- make_track(
  polarBear_GPS,
  .x = Longitude,
  .y = Latitude,
  .t = Datetime,
  id = Bear, 
  crs = 4326)
  
# This is the sf version 
# Linking to GEOS 3.12.2, GDAL 3.9.2, PROJ 9.4.1; sf_use_s2() is

head(trk1)
# A tibble: 6 × 4
     x_    y_ t_                  id    
* <dbl> <dbl> <dttm>              <chr> 
1 -93.2  58.8 2019-08-26 20:20:11 X17517
2 -93.2  58.8 2019-08-26 20:25:10 X17517
3 -93.2  58.8 2019-08-26 20:30:09 X17517
4 -93.2  58.8 2019-08-26 20:35:09 X17517
5 -93.2  58.8 2019-08-26 20:40:09 X17517
6 -93.2  58.8 2019-08-26 20:45:09 X17517`

I am also certain this data is WGS84.


Solution

  • You didn't state how you tried to return the CRS, but I assume it was something like st_crs(trk1)? The result of make_track() does not produce an sf object, so you need to use something like attr() to return the CRS:

    library(sf)
    
    # Return object class
    class(trk1)
    # [1] "track_xyt"  "track_xy"   "tbl_df"     "tbl"        "data.frame"
    
    # Try with st_crs()
    st_crs(trk1)
    # Coordinate Reference System: NA
    
    # Now find the name of the variables of interest and where they are located
    str(trk1)
    # trck_xyt [118,858 × 5] (S3: track_xyt/track_xy/tbl_df/tbl/data.frame)
    #  $ x_  : num [1:118858] -93.2 -93.2 -93.2 -93.2 -93.2 ...
    #  $ y_  : num [1:118858] 58.8 58.8 58.8 58.8 58.8 ...
    #  $ t_  : POSIXct[1:118858], format: "2019-08-26 20:20:11" "2019-08-26 20:25:10" ...
    #  $ id  : chr [1:118858] "X17517" "X17517" "X17517" "X17517" ...
    #  $ tod_: Factor w/ 2 levels "day","night": 2 2 2 2 2 2 2 2 2 2 ...
    #  - attr(*, "crs_")=List of 2
    #   ..$ input: chr "EPSG:4326"
    #   ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic Sys"| __truncated__
    #   ..- attr(*, "class")= chr "crs"
    
    # Return EPSG code
    attr(trk1, "crs_")$input
    # [1] "EPSG:4326"
    
    # Return full WKT 
    attr(trk1, "crs_")$wkt
    # [1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"