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.
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]]"