jmsigner / amt

38 stars 13 forks source link

invalid crs from Move object #39

Closed coverton-usgs closed 2 years ago

coverton-usgs commented 3 years ago

I had a workflow using a single individual move object (originally subset from a larger moveStack).

I was able to use hr_mcp, hr_locoh just fine, but received an error with hr_kde:

bird <- allmovebankdata[["X193540.rail"]]
library(amt)
birdtrack <- as_track(bird)
mcp <- hr_mcp(birdtrack)
locohK <- hr_locoh(birdtrack)
kde2 <- hr_kde(birdtrack, levels = c(0.5, 0.95))

Error in sp::CRS(SRS_string = from$wkt) : 
  no arguments in initialization list

also with

kde2 <- birdtrack %>% hr_kde(h=hr_kde_pi(.),trast = raster(as_sp(.), nrow = 1000, ncol = 1000), levels = c(0.5, 0.95), rand_buffer = 10)
Error in checkSlotAssignment(object, name, value) : 
  assignment of an object of class “character” is not valid for slot ‘proj4string’ in an object of class “Spatial”; is(value, "CRS") is not TRUE

When I make a track without using 'as_track' or assigning a crs it works.

bird2 <- as(bird,"data.frame")
birdtrack2 <- make_track(bird2,.x= location_long, .y = location_lat, .t = timestamp)
kde <- hr_kde(birdtrack2)
class(kde)
[1] "kde"     "hr_prob" "hr"      "list"  

which is great, but my data is in lat-long so hr_area produces incomparable hr sizes.

coverton-usgs commented 3 years ago

Maybe this is an 'sp' issue?

> data(deer)
> hr_kde(deer)
Error in sp::CRS(SRS_string = from$wkt) : 
  no arguments in initialization list
In addition: Warning message:
In CPL_crs_from_input(x) :
  GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.
> sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] amt_0.1.4        move_4.0.4       rgdal_1.4-4      raster_3.4-5     sp_1.4-5         geosphere_1.5-10

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5         compiler_3.6.0     pillar_1.4.6       class_7.3-15       tools_3.6.0        digest_0.6.25      checkmate_1.9.4   
 [8] lubridate_1.7.4    memoise_1.1.0      lifecycle_0.2.0    tibble_3.0.3       lattice_0.20-38    pkgconfig_2.0.3    rlang_0.4.7       
[15] Matrix_1.2-17      bibtex_0.4.2.2     DBI_1.0.0          cli_2.0.2          rstudioapi_0.11    parallel_3.6.0     e1071_1.7-2       
[22] httr_1.4.2         stringr_1.4.0      dplyr_1.0.2        xml2_1.3.2         rgeos_0.4-3        generics_0.0.2     vctrs_0.3.4       
[29] classInt_0.4-3     tidyselect_1.1.0   gbRd_0.4-11        grid_3.6.0         glue_1.4.2         sf_0.9-8           R6_2.4.1          
[36] fansi_0.4.1        Rdpack_0.11-1      survival_3.2-7     tidyr_1.1.2        purrr_0.3.4        magrittr_1.5       units_0.6-3       
[43] backports_1.1.9    codetools_0.2-16   ellipsis_0.3.1     splines_3.6.0      assertthat_0.2.1   utf8_1.1.4         KernSmooth_2.23-15
[50] stringi_1.5.3      crayon_1.3.4   
coverton-usgs commented 3 years ago

Or maybenot - the original track produced using as_track had a valid CRS when used in hr_mcp

> get_crs(mcp)
Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
  wkt:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ID["EPSG",6326]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8901]],
            CS[ellipsoidal,2],
                AXIS["longitude",east,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]],
                AXIS["latitude",north,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
        METHOD["Geocentric translations (geog2D domain)",
            ID["EPSG",9603]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]]]]
jmsigner commented 3 years ago

Could you try the development version of amt, I think I fixed this already. You can install it with remotes::install_github("jmsigner/amt").

coverton-usgs commented 3 years ago

The development version did produce results but I am still not getting exactly what was expected. There appears to be an output unit difference between hr_kde and hr_mcp and hr_locoh - when changing the projection I get almost but not quite the same output hr_area() for the mcp; a changed output unit for hr_locoh(); and an error for hr_kde().

> birdtrack <- as_track(bird)
Warning message:
In sp::proj4string(x) : CRS object has comment, which is lost in output
> get_crs(birdtrack)
Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs 
  wkt:
GEOGCRS["unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
> kde <- hr_kde(birdtrack, levels = 0.95)
> mcp <- hr_mcp(birdtrack, levels = 0.95)
> locohK_10 <- hr_locoh(birdtrack, levels = 0.95)
> 
> hr_area(kde, units=TRUE)
# A tibble: 1 x 3
  level what           area
  <dbl> <fct>         <dbl>
1  0.95 estimate 0.00000323
> hr_area(mcp, units=TRUE)
# A tibble: 1 x 3
  level what         area
  <dbl> <fct>       [m^2]
1  0.95 estimate 42553.36
> hr_area(locohK_10, units=TRUE)
# A tibble: 1 x 3
  level what           area
  <dbl> <fct>         <dbl>
1  0.95 estimate 0.00000136
> 
> birdtrack2 <- transform_coords(birdtrack, crs_to = 32610)
> get_crs(birdtrack2)
[1] 32610
> kde2 <- hr_kde(birdtrack2, levels = 0.95)
Error in as(get_crs(x), "CRS") : 
  no method or default for coercing “numeric” to “CRS”
> mcp2 <- hr_mcp(birdtrack2, levels = 0.95)
> locohK_102 <- hr_locoh(birdtrack2, levels = 0.95)
> 
> hr_area(kde2, units=TRUE)
Error in hr_area(kde2, units = TRUE) : object 'kde2' not found
> hr_area(mcp2, units=TRUE)
# A tibble: 1 x 3
  level what         area
  <dbl> <fct>       [m^2]
1  0.95 estimate 42655.53
> hr_area(locohK_102, units=TRUE)
# A tibble: 1 x 3
  level what       area
  <dbl> <fct>     <dbl>
1  0.95 estimate 11242.

plot(kde) image

plot(mcp) image

plot(locohK_10) image

jmsigner commented 2 years ago

This should all work as expected now:

> library(amt)
> data(deer)
> mcp <- hr_mcp(deer)
> kde <- hr_kde(deer)
> akde <- hr_akde(deer)
> locoh <- hr_locoh(deer)
> hr_area(mcp, units = TRUE)
# A tibble: 1 × 3
  level what          area
  <dbl> <chr>        [m^2]
1  0.95 estimate 30896346.
> hr_area(kde, units = TRUE)
# A tibble: 1 × 3
  level what          area
  <dbl> <chr>        [m^2]
1  0.95 estimate 42029602.
> hr_area(akde, units = TRUE)
# A tibble: 3 × 3
  level what            area
  <dbl> <chr>          [m^2]
1  0.95 lci (0.95) 30263901.
2  0.95 estimate   32816315.
3  0.95 uci (0.95) 34687178.
> hr_area(locoh, units = TRUE)
# A tibble: 1 × 3
  level what         area
  <dbl> <chr>       [m^2]
1  0.95 estimate 3597168.