ateucher / rmapshaper

An R wrapper for the mapshaper javascript library
http://andyteucher.ca/rmapshaper/
Other
200 stars 13 forks source link

Issue with oversimplification of features with ms_innerlines() #161

Closed katherinesiegel closed 1 year ago

katherinesiegel commented 1 year ago

Hi! This package is super helpful -- thanks for your work putting it together! I'm getting an error with ms_innerlines() that I was hoping you could shed some light on. I have a sf dataframe test (attached at the bottom of this post as a csv) and I want to get the internal boundaries between two multipolygons:

bound_test <- test %>% ms_innerlines() %>% as_tibble() %>% st_as_sf()

When I run this, I get the following error message "Error: Cannot convert result to an sf object. It is likely too much simplification was applied and all features were reduced to null." I also tried it on a different set of multipolygons and got the same error. Any idea what that means/how to fix it?

csv here: test_github.csv

Thanks!

ateucher commented 1 year ago

Thanks for the report! That error message is clearly not helpful to you here as you're not using ms_simplify(). That error is thrown when there are no geometries in the result (which almost always results from oversimplification in ms_simplify()).

Can you please provide the example as a geojson file? (sf::st_write(test, "test_github.geojson")). Parsing the geometry in csv column is a bit tricky.

katherinesiegel commented 1 year ago

Thanks for the speedy response! I wasn't able to attach the file as a geojson here, but here's a link to a repo with the file.

ateucher commented 1 year ago

Are you using the latest version of rmapshaper (0.5.0)? I don't get an error, but I do get a result with zero geometries (so the issue is real):

library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(rmapshaper)

poly <- st_read("https://raw.githubusercontent.com/katherinesiegel/ms_innerlines/master/test_github.geojson")
#> Reading layer `test_github' from data source 
#>   `https://raw.githubusercontent.com/katherinesiegel/ms_innerlines/master/test_github.geojson' 
#>   using driver `GeoJSON'
#> Simple feature collection with 2 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -1016202 ymin: 1829426 xmax: -1011419 ymax: 1835521
#> Projected CRS: NAD83 / Conus Albers

ms_innerlines(poly)
#> Simple feature collection with 0 features and 0 fields
#> Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
#> Projected CRS: NAD83 / Conus Albers
#> [1] geometry
#> <0 rows> (or 0-length row.names)

However it is the underlying mapshaper library that is unable to extract the shared boundaries. Using the command-line mapshaper (note the message about shared boundaries):

$ wget https://raw.githubusercontent.com/katherinesiegel/ms_innerlines/master/test_github.geojson
$ mapshaper -i test_github.geojson -innerlines -o foo.geojson         
[innerlines] No shared boundaries were found
[o] RFC 7946 warning: non-WGS84 GeoJSON output.
[o] Wrote foo.geojson
$ cat foo.geojson 
{"type":"GeometryCollection", "geometries": [
]}

It appears that the boundaries between polygons are not quite close enough to be detected as "shared". I haven't implemented ms_snap() to snap boundaries together, but you can use ms_simplify() to do it - set the keep = 1 to keep all vertices (i.e., apply no simplification), but play with the snap_interval parameter to get the appropriate amount of shared boundary snapping. e.g.:

library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(rmapshaper)

poly <- st_read("https://raw.githubusercontent.com/katherinesiegel/ms_innerlines/master/test_github.geojson", 
                quiet = TRUE)

ms_simplify(poly, keep = 1, snap_interval = 10) |> 
ms_innerlines() |> 
  plot()

Created on 2023-06-09 with reprex v2.0.2

ateucher commented 1 year ago

Another note - if you zoom right in in a GIS or with mapview::mapview() you can see the gaps between the boundaries.

image
katherinesiegel commented 1 year ago

Wow, this is so helpful, thank you so much! I see what you mean about the polygons not actually being adjacent...whoops! When I ran the code with ms_simplify() before ms_innerlines(), it worked.

To answer your question about versions, I'm using rmapshaper 0.5.0 and sf 1.0-13 and still getting the error message about too much simplification when I try to run the original code chunk (bound_test <- test %>% ms_innerlines() %>% as_tibble() %>% st_as_sf()). So it's a little puzzling that I'm getting the error message but it doesn't show up on your end!

Maybe it's because some of the dependencies for sf that I'm running are older? I don't know much/anything about the underlying libraries, but here's what I get when I load sf:

`library(sf

> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE`

Anyway, thanks again so much for your help and for putting together this super helpful package!

ateucher commented 1 year ago

Maybe it's because some of the dependencies for sf that I'm running are older

That could be - I need to clean the handling of this kind of thing anyway.

Happy to help!