joelgombin / concaveman

A very fast 2D concave hull algorithm
https://joelgombin.github.io/concaveman/
66 stars 7 forks source link

Possible error in algorithm - polygon fails to include all points #17

Open AdamCSmithCWS opened 1 year ago

AdamCSmithCWS commented 1 year ago

Thanks so much for your awesome package. I've been using it a lot. I've been perplexed lately by some odd situations where the algorithm seems to ignore one or a few points when it's creating the polygon.

I've uploaded a small dataset of coordinates to the issue and I hope the code below provides a full reproducible example.

coordinates_example.csv

I've seen this issue with a number of point datasets recently. The data represent survey locations from a continental scale bird monitoring program, but for this example the geographic aspects of the points are irrelevant (I think).

library(tidyverse)
library(sf)
library(concaveman)
# import text file of coordinates
coords <- read.csv("coordinates_example.csv") %>% 
  st_as_sf(coords = c("X","Y"))

Running the concave() function on these points creates a polygon that does not intersect all points.

concave_poly <- concaveman(coords) %>% 
  mutate(in_poly = TRUE) #adds a column for following join

missed <- coords %>% 
  st_join(.,concave_poly,
          join = st_intersects) %>% 
  filter(is.na(in_poly))

In this case, there are three points that are outside of the polygon. However, two of them are extremely close, such that a 1-unit buffer on the polygon includes them. So not a problem in this context, likely more an issue with the precision of the st_intersects() function. So I added a 1-unit buffer below, but there is still one point outside the polygon, and it's far outside. It seems to be completely ignored by the algorithm.


concave_poly <- concaveman(coords) %>% 
  st_buffer(.,1) %>% #just adds a 1-unit buffer - simplifies st_join
  mutate(in_poly = TRUE) #adds a column for following join

missed <- coords %>% 
  st_join(.,concave_poly,
          join = st_intersects) %>% 
  filter(is.na(in_poly))

# plot the polygon with the missing point in red
display <- ggplot()+
  geom_sf(data = concave_poly)+ #concave polygon
  geom_sf(data = coords,alpha = 0.3)+ #original points
  geom_sf(data = missed,colour = "red")+ #points that do not intersect polygon
  geom_sf_text(data = missed,
               aes(label = point),
               size = 3)+
  labs(title = "red point is far outside the polygon")

display

This plot shows the red point that's completely missed by the polygon. plot_polygon_overlap

sessionInfo()
R version 4.2.0 (2022-04-22 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

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

other attached packages:
 [1] concaveman_1.1.0 sf_1.0-9         forcats_0.5.1   
 [4] stringr_1.5.0    dplyr_1.1.0      purrr_0.3.4     
 [7] readr_2.1.2      tidyr_1.2.0      tibble_3.1.7    
[10] ggplot2_3.4.0    tidyverse_1.3.1 

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3       lubridate_1.8.0    class_7.3-20      
 [4] assertthat_0.2.1   digest_0.6.29      utf8_1.2.2        
 [7] V8_4.2.2           R6_2.5.1           cellranger_1.1.0  
[10] backports_1.4.1    reprex_2.0.1       evaluate_0.15     
[13] e1071_1.7-11       httr_1.4.3         pillar_1.7.0      
[16] rlang_1.0.6        curl_4.3.2         readxl_1.4.0      
[19] rstudioapi_0.13    rmarkdown_2.14     munsell_0.5.0     
[22] proxy_0.4-27       broom_0.8.0        compiler_4.2.0    
[25] modelr_0.1.8       xfun_0.31          pkgconfig_2.0.3   
[28] htmltools_0.5.2    tidyselect_1.2.0   fansi_1.0.3       
[31] crayon_1.5.1       tzdb_0.3.0         dbplyr_2.2.1      
[34] withr_2.5.0        grid_4.2.0         jsonlite_1.8.0    
[37] gtable_0.3.0       lifecycle_1.0.3    DBI_1.1.3         
[40] magrittr_2.0.3     units_0.8-0        scales_1.2.0      
[43] KernSmooth_2.23-20 cli_3.6.0          stringi_1.7.6     
[46] farver_2.1.0       fs_1.5.2           xml2_1.3.3        
[49] ellipsis_0.3.2     generics_0.1.2     vctrs_0.5.2       
[52] tools_4.2.0        glue_1.6.2         hms_1.1.1         
[55] fastmap_1.1.0      yaml_2.3.5         colorspace_2.0-3  
[58] classInt_0.4-7     rvest_1.0.2        knitr_1.39        
[61] haven_2.5.0 
AdamCSmithCWS commented 1 year ago

FYI, I updated everything to R 4.3, and the same error occurs:

sessionInfo()
R version 4.3.0 (2023-04-21 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

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

other attached packages:
 [1] concaveman_1.1.0 sf_1.0-12        lubridate_1.9.2  forcats_1.0.0    stringr_1.5.0   
 [6] dplyr_1.1.2      purrr_1.0.1      readr_2.1.4      tidyr_1.3.0      tibble_3.2.1    
[11] ggplot2_3.4.2    tidyverse_2.0.0 

loaded via a namespace (and not attached):
 [1] utf8_1.2.3         generics_0.1.3     class_7.3-21       KernSmooth_2.23-20
 [5] stringi_1.7.12     hms_1.1.3          digest_0.6.31      magrittr_2.0.3    
 [9] evaluate_0.20      grid_4.3.0         timechange_0.2.0   fastmap_1.1.1     
[13] jsonlite_1.8.4     e1071_1.7-13       DBI_1.1.3          fansi_1.0.4       
[17] scales_1.2.1       cli_3.6.1          rlang_1.1.0        units_0.8-1       
[21] munsell_0.5.0      withr_2.5.0        yaml_2.3.7         tools_4.3.0       
[25] tzdb_0.3.0         colorspace_2.1-0   curl_5.0.0         vctrs_0.6.2       
[29] R6_2.5.1           proxy_0.4-27       lifecycle_1.0.3    classInt_0.4-9    
[33] V8_4.3.0           pkgconfig_2.0.3    pillar_1.9.0       gtable_0.3.3      
[37] glue_1.6.2         Rcpp_1.0.10        xfun_0.39          tidyselect_1.2.0  
[41] rstudioapi_0.14    knitr_1.42         farver_2.1.1       htmltools_0.5.5   
[45] rmarkdown_2.21     compiler_4.3.0   
hunterdanielsmith commented 1 year ago

I'm seeing a similar issue.

bad_boundary

cysouw commented 5 months ago

I see the same error. I did some testing to see what influences the omission of points. Re-ordering of the points did not change the omission. Adding some jitter to the values did not solve it either.

However, when I manually round the coordinates of the problematic point just slightly (in my case rounding the longitude from 8.34688 to 8.35), then the point is included in the hull. Note that rounding to 8.347 is not enough to get the point included.

very strange bug...