search

jmlondon / pathroutr

R package for routing tracks around barrier polygons
https://jmlondon.github.io/pathroutr/
16 stars 2 forks source link

Error: arguments have different crs #7

Closed jamesgrecian closed 1 year ago

jamesgrecian commented 3 years ago

When using the pathroutr package to pull some simulated tracks out to sea I often get the following: Error: arguments have different crs I can't work out the origin of the error as the land polygon and movement path have the same crs info. Any help would be greatly appreciated!

Here's a reprex that should recreate the error using a small sample of my data, sessionInfo below.

# load libraries
require(foieGras)
#> Loading required package: foieGras
#> Registered S3 method overwritten by 'spatstat.geom':
#>   method     from
#>   print.boxx cli
require(pathroutr)
#> Loading required package: pathroutr
require(sf)
#> Loading required package: sf
#> Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
require(tidyverse)
#> Loading required package: tidyverse
require(rnaturalearth)
#> Loading required package: rnaturalearth

# load world shapefile for visgraph generation
world_mc <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  st_transform(., crs = "+proj=merc +units=km +datum=WGS84 +no_defs") %>%
  sf::st_make_valid()

# example simulated track
trs_df <- structure(list(id = c("168556", "168556", "168556", "168556", "168556",
                                "168556", "168556", "168556", "168556", "168556", 
                                "168556", "168556", "168556", "168556", "168556",
                                "168556", "168556", "168556", "168556", "168556",
                                "168556", "168556", "168556", "168556", "168556"),
                         model = c(`168556` = "rw", `168556` = "rw", `168556` = "rw", `168556` = "rw", 
                                   `168556` = "rw", `168556` = "rw", `168556` = "rw", `168556` = "rw", 
                                   `168556` = "rw", `168556` = "rw", `168556` = "rw", `168556` = "rw", 
                                   `168556` = "rw", `168556` = "rw", `168556` = "rw", `168556` = "rw", 
                                   `168556` = "rw", `168556` = "rw", `168556` = "rw", `168556` = "rw", 
                                   `168556` = "rw", `168556` = "rw", `168556` = "rw", `168556` = "rw", 
                                   `168556` = "rw"),
                         rep = c(12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
                                 12, 12, 12, 12, 12, 12, 12, 12, 12, 12),
                         date = structure(c(1553515200, 1553601600, 1553688000, 1553774400, 1553860800,
                                            1553947200, 1554033600, 1554120000, 1554206400, 1554292800, 
                                            1554379200, 1554465600, 1554552000, 1554638400, 1554724800, 
                                            1554811200, 1554897600, 1554984000, 1555070400, 1555156800,
                                            1555243200, 1555329600, 1555416000, 1555502400, 1555588800),
                                          tzone = "", class = c("POSIXct", "POSIXt")),
                         lon = c(-58.5659198518063, -59.4563814383551, -60.6391495534974, -60.7999008399566, -61.4353755608006,
                                 -60.1036452713862, -60.7637828117047, -60.1314427063631, -61.4144372028938, -62.3209388971776,
                                 -62.5137719397771, -63.7938729728152, -63.3769100131878, -66.3788977685044, -70.167473931186,
                                 -69.096705957986, -71.6422524443119, -70.7455427679063, -69.7022629494988, -71.5636192570846,
                                 -70.0303993492376, -70.2205830348059, -67.2596659138549, -63.2979622571531, -62.3932336178804),
                         lat = c(48.54721903525, 48.2493851506603, 48.2798033736316, 49.1932412435527, 50.2664408417163,
                                 49.1331710713557, 49.9923879043313, 49.4527006377747, 49.6450143733593, 48.8579857940244,
                                 48.6432001089196, 47.5800419966615, 47.4319873554999, 47.7755617702211, 48.0893231015444,
                                 47.6499005537445, 48.1651413299542, 47.4948035120105, 47.0469874436098, 47.7349508799141, 
                                 46.6985786087455, 48.0652674724875, 47.5644464599778, 47.6271699870094, 46.8801193363709),
                         x = c(-6519.52837574275, -6618.65410612834, -6750.3192504325, -6768.21400178549, -6838.95472412184,
                               -6690.70718643026, -6764.19336127204, -6693.80158273724, -6836.62387678161, -6937.53518379253,
                               -6959.00125990282, -7101.50145506457, -7055.08535071919, -7389.26509900867, -7811.00746826993,
                               -7691.81012273556, -7975.17906138397, -7875.35779681709, -7759.22041867707, -7966.42565502237,
                               -7795.74839560672, -7816.91954664138, -7487.31176045603, -7046.29692671814, -6945.5829952882),
                         y = c(6166.33807868558, 6116.55011251715, 6121.62161758653, 6275.35417629948, 6459.66817858585,
                               6265.15735729721, 6412.21007861329, 6319.54068410451, 6352.44414098663, 6218.60322059442,
                               6182.44568835735, 6005.71229363521, 5981.38876931086, 6037.94002922736, 6089.91348957851,
                               6017.21321130505, 6102.52030690898, 5991.70024875042, 5918.45858575017, 6031.23608432989,
                               5861.90218570855, 6085.9175029917, 6003.1468861643, 6013.46934743584, 5891.32516405129)), 
                    row.names = c(NA, -25L),
                    class = c("tbl_df", "tbl", "data.frame"))

# convert track to sf for pathroutr
trs_sf <- trs_df %>%
  sf::st_as_sf(coords = c("x", "y"), crs = "+proj=merc +units=km +datum=WGS84 +no_defs")

# generate land region
land_region <- sf::st_buffer(trs_sf, dist = 50) %>% 
  sf::st_union() %>% 
  sf::st_convex_hull() %>% 
  sf::st_intersection(world_mc) %>% 
  st_collection_extract('POLYGON') %>% 
  st_sf()

# create visgraph
vis_graph <- prt_visgraph(land_region)

# pulling points to land should fail for at least one of the simulated tracks...
track_pts_fix <- trs_sf %>%
  prt_trim(land_region) %>%
  prt_reroute(., land_region, vis_graph)
#> Error: arguments have different crs

Created on 2021-05-20 by the reprex package (v2.0.0)

> traceback()
8: stop("arguments have different crs", call. = FALSE)
7: rbind(deparse.level, ...)
6: rbind(degenerate_paths, regular_paths)
5: stplanr::sum_network_routes(sln = vis_graph, start = segs_tbl$start_node, 
       end = segs_tbl$end_node)
4: prt_shortpath(., vis_graph)
3: get_barrier_segments(trkpts, barrier) %>% prt_nearestnode(vis_graph) %>% 
       prt_shortpath(vis_graph)
2: prt_reroute(., land_region, vis_graph)
1: trs_sf %>% prt_trim(land_region) %>% prt_reroute(., land_region, 
       vis_graph)
> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

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

other attached packages:
 [1] rnaturalearth_0.1.0 sf_0.9-8            pathroutr_0.1.1     foieGras_0.7-6.9267 forcats_0.5.1      
 [6] stringr_1.4.0       dplyr_1.0.6         purrr_0.3.4         readr_1.4.0         tidyr_1.1.3        
[11] tibble_3.1.2        ggplot2_3.3.3       tidyverse_1.3.1     reprex_2.0.0       

loaded via a namespace (and not attached):
  [1] colorspace_2.0-1        deldir_0.2-10           ellipsis_0.3.2          class_7.3-18           
  [5] trip_1.8.5              fs_1.5.0                rstudioapi_0.13         proxy_0.4-25           
  [9] spatstat.data_2.1-0     fansi_0.4.2             mvtnorm_1.1-1           lubridate_1.7.10       
 [13] xml2_1.3.2              codetools_0.2-18        splines_4.0.3           knitr_1.33             
 [17] polyclip_1.10-0         jsonlite_1.7.2          reproj_0.4.2            broom_0.7.6            
 [21] dbplyr_2.1.1            rgeos_0.5-5             spatstat.sparse_2.0-0   clipr_0.7.1            
 [25] compiler_4.0.3          httr_1.4.2              backports_1.2.1         assertthat_0.2.1       
 [29] Matrix_1.3-2            gmm_1.6-6               cli_2.5.0               htmltools_0.5.1.1      
 [33] tools_4.0.3             igraph_1.2.6            rnaturalearthdata_0.1.0 gtable_0.3.0           
 [37] glue_1.4.2              Rcpp_1.0.6              cellranger_1.1.0        raster_3.4-10          
 [41] vctrs_0.3.8             nlme_3.1-152            tmvtnorm_1.4-10         xfun_0.22              
 [45] ps_1.6.0                rvest_1.0.0             lifecycle_1.0.0         goftest_1.2-2          
 [49] crsmeta_0.3.0           MASS_7.3-53.1           zoo_1.8-9               scales_1.1.1           
 [53] spatstat.core_2.1-2     hms_1.1.0               spatstat.utils_2.1-0    parallel_4.0.3         
 [57] proj4_1.0-10.1          sandwich_3.0-0          TMB_1.7.20              curl_4.3.1             
 [61] yaml_2.2.1              geosphere_1.5-10        rpart_4.1-15            traipse_0.2.0          
 [65] stringi_1.6.2           highr_0.9               maptools_1.1-1          e1071_1.7-6            
 [69] boot_1.3-27             rlang_0.4.11            pkgconfig_2.0.3         evaluate_0.14          
 [73] lattice_0.20-41         tensor_1.5              patchwork_1.1.1         processx_3.5.1         
 [77] tidyselect_1.1.1        magrittr_2.0.1          R6_2.5.0                generics_0.1.0         
 [81] DBI_1.1.1               foreign_0.8-81          pillar_1.6.1            haven_2.4.1            
 [85] withr_2.4.2             mgcv_1.8-35             units_0.7-1             abind_1.4-5            
 [89] sp_1.4-5                modelr_0.1.8            crayon_1.4.1            KernSmooth_2.23-18     
 [93] utf8_1.2.1              spatstat.geom_2.1-0     nabor_0.5.0             rmarkdown_2.7          
 [97] grid_4.0.3              readxl_1.3.1            callr_3.7.0             CircStats_0.2-6        
[101] digest_0.6.27           classInt_0.4-3          stats4_4.0.3            munsell_0.5.0          
[105] stplanr_0.8.2
jmlondon commented 3 years ago

James

Thank you for filing this issue and, especially, for the reprex.

I think I have an answer for what is happening, but I don't fully understand why.

The issue appears tied to your use of a custom proj4string (+proj=merc +units=km +datum=WGS84 +no_defs). I've seen recent discussions regarding changes in sf and other spatial packages that, now, support newer versions of PROJ. My understanding was that proj4strings were discouraged for describing crs but would still be supported. Maybe your use of units=km is no longer supported?? I'll try to so some additional reading and sleuthing to learn more on this topic. There should be a way to specify a custom proj4string as you do.

Regardless, though, the following adjustments to your code seem to work. I'm simply relying on the ESPG code of 3857 for mercator instead of your proj4string. Note, this creates geometry that is the same as your x/y coordinates but in meters.

# load world shapefile for visgraph generation
world_mc <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  st_transform(., crs = 3857) %>%
  sf::st_make_valid()
# convert track to sf for pathroutr
trs_sf <- trs_df %>%
  sf::st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  sf::st_transform(3857)
jamesgrecian commented 3 years ago

Thanks @jmlondon that's solved the issue, although as you say it's a strange error given the crs info is the same.

makratofil commented 3 years ago

@jmlondon I know this issue was closed, but if it helps at all, I ran into the same issue specifying "crs = 102007" (Hawaii Albers Equal Area Conic). I've used this projection before when running pathroutr but didn't run into any issues regarding crs' until today. Using the ESPG 3857 code solved the problem.

R version 4.0.5 (2021-03-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

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] stplanr_0.8.2    pathroutr_0.1.1  concaveman_1.1.0 janitor_2.1.0    ggspatial_1.1.5  ggplot2_3.3.3    crawl_2.2.3      mapview_2.9.9   
 [9] sf_0.9-8         purrr_0.3.4      lubridate_1.7.10 dplyr_1.0.6     

loaded via a namespace (and not attached):
 [1] tidyr_1.1.3             jsonlite_1.7.2          shiny_1.6.0             assertthat_0.2.1        sp_1.4-5                stats4_4.0.5           
 [7] yaml_2.2.1              pillar_1.6.1            lattice_0.20-41         glue_1.4.2              uuid_0.1-4              digest_0.6.27          
[13] promises_1.2.0.1        snakecase_0.11.0        colorspace_2.0-1        leaflet.providers_1.9.0 htmltools_0.5.1.1       httpuv_1.6.1           
[19] pkgconfig_2.0.3         raster_3.4-10           xtable_1.8-4            mvtnorm_1.1-1           scales_1.1.1            webshot_0.5.2          
[25] brew_1.0-6              svglite_2.0.0           satellite_1.0.2         later_1.2.0             tibble_3.1.1            proxy_0.4-25           
[31] farver_2.1.0            generics_0.1.0          ellipsis_0.3.2          withr_2.4.2             cli_2.5.0               magrittr_2.0.1         
[37] crayon_1.4.1            mime_0.10               maptools_1.1-1          fansi_0.4.2             foreign_0.8-81          class_7.3-18           
[43] tools_4.0.5             geosphere_1.5-10        nabor_0.5.0             lifecycle_1.0.0         stringr_1.4.0           V8_3.4.2               
[49] munsell_0.5.0           packrat_0.5.0           compiler_4.0.5          e1071_1.7-6             systemfonts_1.0.2       tinytex_0.31           
[55] rlang_0.4.11            classInt_0.4-3          units_0.7-1             grid_4.0.5              leafpop_0.1.0           rstudioapi_0.13        
[61] htmlwidgets_1.5.3       igraph_1.2.6            crosstalk_1.1.1         proj4_1.0-10.1          leafem_0.1.6            base64enc_0.1-3        
[67] gtable_0.3.0            codetools_0.2-18        curl_4.3.1              DBI_1.1.1               R6_2.5.0                rgeos_0.5-5            
[73] fastmap_1.1.0           utf8_1.2.1              KernSmooth_2.23-18      stringi_1.5.3           Rcpp_1.0.6              vctrs_0.3.8            
[79] png_0.1-7               leaflet_2.0.4.1         tidyselect_1.1.1        xfun_0.22
jmlondon commented 3 years ago

Thanks for reporting this, @makratofil ... I'm going to re-open this issue and create a new issue to track progress.

jmlondon commented 3 years ago

@makratofil one thing to try, in your case, is to specify the CRS as 'ESRI:102007' and see if that works.

makratofil commented 3 years ago

@jmlondon just got back to this analysis, and the 'ESRI:102007' works! thanks!

jmlondon commented 3 years ago

@makratofil and @jamesgrecian ... I just wanted to follow up and let you know that I think the recent updates I've made (now, in version 0.2.1) may have fixed these issues.

I still don't understand why these issues were happening, but I suspect the switch to {sfnetworks} from {stplanr} may have had the added benefit of fixing this issue.

If either of you have a chance to test and confirm on your end, then I'll close this issue.