Polygons disappear after ms_erase

[sandra-ab]

Hi,

First of all thanks for rmapshaper, it is a wonderfully handy package.

I have been occasionally encountering a quirk using ms_erase and can’t really figure out what is happening.

For context, I am using the function in a habitat mapping toolkit, and users can supply their own spatial file of hedgerows that then get added to a vector land cover map (OS MasterMap). The (simplified) workflow for this is:

1. Buffer the linear features (hedgerows) to turn them into polygons (they get unioned and validated to be as tidy as possible, even though I have no control on what the user supplies…)
2. Erase (using ms_erase) the corresponding area from the basemap polygons
3. Combine the two layers

Now this works beautifully almost every time (> 99.9%), but once in a while the returned output of ms_erase is missing part of a polygon, leaving a gaping hole in my map. In the example here, I also noticed that this does not happen when I use the system rmapshaper instead. The plot shows my two input layers (left), the problem (center), and the expected correct outcome (right; using system rmapshaper or QGIS).

Sorry I couldn’t reproduce the problem with simpler data, but as I said it only happens very rarely - and I can't figure out what combination of features sets it off (all input polygons are valid sf single-part polygons). I am attaching the two layers shown in the plots: hedgemap.zip

Any pointers on what the problem might be would be greatly appreciated. I would prefer not forcing my users to install Node.js.

My code:

library(rmapshaper)
library(sf)

# Import data
unzip("hedgemap.zip")
hedge <- st_read("hedgemap/hedge.geojson")
poly <- st_read("hedgemap/poly.geojson")

# Erase hedgerows from polygons
wrong <- ms_erase(poly, hedge)
right <- ms_erase(poly, hedge, sys = TRUE)

And the output of sessionInfo():

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=French_Canada.1252  LC_CTYPE=French_Canada.1252    LC_MONETARY=French_Canada.1252
[4] LC_NUMERIC=C                   LC_TIME=French_Canada.1252    

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

other attached packages:
[1] sf_0.9-8         rmapshaper_0.4.4

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.6         lattice_0.20-41    tidyr_1.1.3        class_7.3-18       digest_0.6.27     
 [6] packrat_0.6.0      utf8_1.2.1         V8_3.4.2           R6_2.5.0           cellranger_1.1.0  
[11] backports_1.2.1    evaluate_0.14      e1071_1.7-6        geojson_0.3.4      ggplot2_3.3.3     
[16] pillar_1.6.1       rlang_0.4.11       lazyeval_0.2.2     curl_4.3.1         readxl_1.3.1      
[21] geojsonlint_0.4.0  data.table_1.14.0  car_3.0-10         geojsonio_0.9.4    rmarkdown_2.8     
[26] jqr_1.2.1          foreign_0.8-81     munsell_0.5.0      proxy_0.4-25       broom_0.7.6       
[31] compiler_4.0.5     xfun_0.22          pkgconfig_2.0.3    rgeos_0.5-5        htmltools_0.5.1.1 
[36] tidyselect_1.1.1   tibble_3.1.1       httpcode_0.3.0     rio_0.5.26         jsonvalidate_1.1.0
[41] fansi_0.4.2        crayon_1.4.1       dplyr_1.0.6        ggpubr_0.4.0       crul_1.1.0        
[46] grid_4.0.5         jsonlite_1.7.2     gtable_0.3.0       lifecycle_1.0.0    DBI_1.1.1         
[51] magrittr_2.0.1     units_0.7-1        scales_1.1.1       KernSmooth_2.23-18 zip_2.1.1         
[56] stringi_1.5.3      carData_3.0-4      ggsignif_0.6.1     sp_1.4-5           ellipsis_0.3.2    
[61] generics_0.1.0     vctrs_0.3.8        openxlsx_4.2.3     geojsonsf_2.0.1    tools_4.0.5       
[66] forcats_0.5.1      glue_1.4.2         purrr_0.3.4        hms_1.0.0          abind_1.4-5       
[71] yaml_2.2.1         colorspace_2.0-1   maptools_1.1-1     rstatix_0.7.0      classInt_0.4-3    
[76] knitr_1.33         haven_2.4.1       

Thank you,
Sandra

[ateucher]

Hi @sandra-ab - thanks for the excellent bug report! I'll investigate shortly

[ateucher]

Hi @sandra-ab - it seems to have been a bug with the bundled version of mapshaper I was using. I've upgraded it, can you test with the development version? You can install with:

# install.packages("remotes")
remotes::install_github("ateucher/rmapshaper", ref = "fix-erase-bug")

[sandra-ab]

Thank you so much @ateucher - works perfectly now!

[ateucher]

That's excellent - thank you for the report!