deleting holes with editor leafpm returns invalid sf

[tim-salabim]

If we use the new leafpm editor to delete holes in polygons, the returned feature is invalid:

library(sf)
library(mapview)
library(leaflet)
library(lwgeom)
# devtools::install_github("r-spatial/mapedit@leafpm")
library(mapedit)
# devtools::install_github("r-spatial/leafpm")
library(leafpm)

outer1 = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
outer2 = matrix(c(11,0,11,1,12,1,12,0,11,0),ncol=2, byrow=TRUE)

pts1 = list(outer1, hole1, hole2)
pts2 = list(outer2)

pl1 = st_sf(geom = st_sfc(st_polygon(pts1)))
pl2 = st_sf(geom = st_sfc(st_polygon(pts2)))

mpl = st_sf(geom = st_combine(rbind(pl1, pl2)))

tst = editFeatures(mpl, editor = "leafpm") # delete one of the holes!

st_is_valid(tst)
# FALSE

tst2 = st_make_valid(tst)
st_is_valid(tst2)
#TRUE

I think we could simply wrap whatever comes back from the editing session in a lwgeom::st_make_valid inside a try(Catch) call as it is not possible to know what users will draw/delete/edit anyway... Thoughts?

[timelyportfolio]

@tim-salabim great idea. Please let me know if I got close to your expectation. Thanks!

[lbusett]

I noticed that this only appears to happen on multipolygons. Casting to POLYGON beforehand and then removing a hole returns valid geometries:

library(mapedit)
library(leafpm)
library(mapedit)
library(sf)

outer1 = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
outer2 = matrix(c(11,0,11,1,12,1,12,0,11,0),ncol=2, byrow=TRUE)

pts1 = list(outer1, hole1, hole2)
pts2 = list(outer2)

pl1 = st_sf(geom = st_sfc(st_polygon(pts1)))
pl2 = st_sf(geom = st_sfc(st_polygon(pts2)))

mpl = st_sf(geom = st_combine(rbind(pl1, pl2)))
tst = editFeatures(mpl, editor = "leafpm") #remove hole
st_is_valid(tst)
>FALSE

mpl2 <- sf::st_cast(mpl, "POLYGON")
tst2 = editFeatures(mpl2, editor = "leafpm") #remove hole 
st_is_valid(tst2)
>TRUE TRUE

Looking at the WKTs returned after removing the hole on the multipolygon:

> st_astext(mpl)
[1] "MULTIPOLYGON(((0 0,10 0,10 10,0 10,0 0),(1 1,1 2,2 2,2 1,1 1),(5 5,5 6,6 6,6 5,5 5)),((11 0,11 1,12 1,12 0,11 0)))"

> st_astext(tst)
[1] "MULTIPOLYGON(((0 0,10 0,10 10,0 10,0 0),(1 1,1 2,2 2,2 1,1 1),(11 0,11 1,12 1,12 0,11 0)))"

it seems that (if I understand WKT properly), we are missing a couple of parenthesis (After "...1, 1 1)" and before "(11 0, 11..", so that the second polygon "appears" to be a hole of the first one. Could this be a bug in the leafpm plugin?

PS: Id'also be careful if possible with st_make_valid: in my experience its results can be difficult to predict.

HTH

[tim-salabim]

Thanks @lbusett! I think you are right that the polygon structure gets mangled up. Indeed seems like an issue in leaflet.pm.

Can you please elaborate on what you mean by

in my experience its results can be difficult to predict.

Any examples where it does something unexpected? I have never experienced any problems other than that some geometries cannot be repaired iirc...

[lbusett]

Hi @tim-salabim. I agree that the "main" problem is in failing to repair geometries. That can be easily "intercepted" (I think @timelyportfolio implemented that, but I wonder if it wouldn't be better to issue an Error, in that case, instead of returning an invalid geometry).

I however think that even in case of success, there are cases where results of st_make_valid should be handled with care. I'll try to explain. For example:

library(lwgeom)
#> Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
outer = matrix(c(0,0,10,0,10,10,11,10, 10,10, 0,10,0,0),ncol=2, byrow=TRUE)
pts   = list(outer)
pl    = st_sf(a = 100, geom = st_sfc(st_polygon(pts)))
plot(pl)

st_is_valid(pl)
#> [1] FALSE
pl2 <- st_make_valid(pl)
st_is_valid(pl2)
#> [1] TRUE
st_astext(pl2)
#> [1] "GEOMETRYCOLLECTION(POLYGON((10 10,10 0,0 0,0 10,10 10)),LINESTRING(10 10,11 10))"

here, in order to repair the geometry we go from a POLYGON to a GEOMETRYCOLLECTION, and a LINESTRING is introduced to maintain the number of vertexes.

library(lwgeom)
#> Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
outer = matrix(c(0,0,10,0,10,10,0,10,10.1,5, 0,0),ncol=2, byrow=TRUE)
pts   = list(outer)
pl    = st_sf(a = 100, geom = st_sfc(st_polygon(pts)))
plot(pl)

pl2 <- st_make_valid(pl)
st_is_valid(pl2)
#> [1] TRUE
st_astext(pl2)
#> [1] "MULTIPOLYGON(((10 5.049505,0 10,10 10,10 5.049505)),((10 4.950495,10 0,0 0,10 4.950495)),((10 4.950495,10 5.049505,10.1 5,10 4.950495)))"

here we go from POLYGON to MULTIPOLYGON, introducing a small "triangle" on the right.

Similar to this second example, but more related to mapedit, it seems to me that at the moment nothing prevents a user from doing something like this:

outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
pts   = list(outer)
pl    = st_sf(a = 100, geom = st_sfc(st_polygon(pts)))
tst = editFeatures(mpl, editor = "leafpm")

Clicking "Done" after editing like this currently returns an invalid geometry, and this can be checked afterwards. However, if st_make_valid is used automatically (as I seem to understand form @timelyportfolio last commit), the edited object would be converted to a valid MULTIPOLYGON without issuing any warnings. The user error would be "masked out", and the invalid geometry transformed into a valid (but "nonsensical") one.

plot(st_make_valid(tst))

Hope this is more or less clear and it makes sense. Apologies if I'm missing anything obvious (or you think these are just "made up" problems unlikely to occur ;-) ).

[timelyportfolio]

@lbusett @tim-salabim thanks so much for the discussion and feedback. Let me see what I can do on the JavaScript side to handle these cases. I don't think there is a perfect solution, but it seems we can do better.

[lbusett]

@timelyportfolio Looking at the demo of leaf.pm (https://geoman.io/studio) it seems that it is possible to set the editor so to prevent self-intersection when creating/editing layers (by setting the allowSelfIntersection "argument" to FALSE) . I think that would be useful (Do not know what would happen however if "loading" a polygon with self-intersected features).

[tim-salabim]

Thanks @lbusett ! Those are great examples that highlight why we should indeed be careful with st_make_valid (especially the second example). My hope is that, as you just alerted to allowSelfIntersection: false will indeed address some of the issues that may happen.

@timelyportfolio in light of these examples, I think we should not force valid geometries, but rather issue a warning if st_is_valid of the edited/created features returns FALSE. This, however, means that we still need to find a suitable way of avoiding invalid geometries when deleting holes from MULTIPOLYGONs.

I think it is acceptable to burden the user with how to deal with invalid geometries. It is a pitty that we don't have a way of getting more information about what it is that makes a geometry invalid. Maybe @edzer knows how to tap into this information. It seems that it should be possible to get more info on what is wrong, as st_intersection for example prints further information along the lines of topology exception, found self-intersection at or near vertex ... Any ideas @edzer?

[edzer]

Something along these lines?

> ring = rbind(c(0,0), c(1,1), c(0,1), c(0,0))
> library(sf)
Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 5.2.0
> st_is_valid(st_polygon(list(ring, ring)))
[1] FALSE
> st_is_valid(st_polygon(list(ring, ring)), reason = TRUE)
[1] "Self-intersection[0 0]"

[tim-salabim]

@edzer perfect! Thanks! I did not know about reason = TRUE

So, it seems to me that it would be better to check whether edited/drawn features are valid and if not, raise a warning with the reason. Thoughts?

[lbusett]

@tim-salabim Makes sense to me.

On a sidenote, I was thinking that it could maybe make sense to have a "repair_geoms" function in mapedit that, given an invalid geometry, prompts the user for interactive correction. Something on this lines:

repair_geoms <- function(sf_object) {
  invalid <- sf_object[!sf::st_is_valid(sf_object), ]
  valid   <- sf_object[sf::st_is_valid(sf_object), ]
  out = editFeatures(invalid, editor = "leafpm")
  if (all(st_is_valid(out))) {
    rbind(valid, invalid)
  } else {
    message("still invalid")
  }
}

what do you think?

[tim-salabim]

Great idea! Fits well into the scope of mapedit. @timelyportfolio what do you think?

[tim-salabim]

@timelyportfolio along the lines of creating valid geometries, can we set style options so that when e.g. a self-intersection occurs the color of the drawn/edited geometry it changes color as in the examples shown here?

[timelyportfolio]

@lbusett @tim-salabim I reverted the commits to check and try to repair validity. Should we save this for the next release? I definitely want to get this solved but worried that it might take some time to get it right. Will wait to hear.

[tim-salabim]

@timelyportfolio i agree that this may be too involved for the next release, especially given the deadline end of Jan. A test whether what has been done produced a valid geometry (together with a warning if it didn't) and having allowSelfIntersection set to false for polygons (if incoming feature is (multi)polygon) should suffice for now.

[timelyportfolio]

@tim-salabim ok will test for valid and warn. I don't see where we can limit allowSelfIntersection to only polygons in Leaflet.pm. For now, do you think it is better to default to FALSE even though it still applies to line or do default to TRUE? Thanks!!

[timelyportfolio]

@tim-salabim @lbusett I upgraded to 2.0.2, and @codeofsumit has made some really nice changes. If you get a chance to test, I'd really appreciate it.

[tim-salabim]

Can't we set the options of leafpm differently depending on what is coming in?

[timelyportfolio]

Yes, but in the case of mixed feature types, we will have to choose one or the other. In the case of mixed, do you prefer FALSE or TRUE.

[tim-salabim]

In case of mixed features, there is one other layer of tests that we could add. If !any(st_diemension(x) == 2) we could also set allowSelfIntersection to true as this means the collection only contains points and lines. In case that returns false, I'd set it to false and issue a warning regarding the selfintersection

[timelyportfolio]

@tim-salabim I tried to add in https://github.com/r-spatial/mapedit/commit/8521dd12b8aa59e8346883de300d79258f01bc75 but it appears the allowSelfIntersection=false does not apply to existing features. I'll do some research to figure this out.

[tim-salabim]

If I do editFeatures(Franconia[5, ], editor = "leafpm") I am not allowed to drag any vertex to create a self intersection which is correct. I even not allowed to delete a vertex - hole or outer ring - that would result in a self intersection. So all seems to work as expected IMO

[lbusett]

Looking good to me! Also the other additions ("moving " button and icons improvements) are quite nice.

Small caveat: on the usual "test polygons with holes" used above, the user can still "eat a hole" doing this:

The resulting geometry is invalid, with the "eaten hole" transformed in a new poly:

I think however that the warning about invalid geoms is sufficient to deal with a case like this: if the user wants to make such a mess, it's his fault ;-)

[tim-salabim]

Took me a while to figure out what you were do there, but that just goes to show that we cannot cover all grounds and need users to behave sensibly. So warnings and errors in the correct places are more important than ever

[lbusett]

Coming back to the original issue, here

https://github.com/lbusett/mapedit/commit/6f127d2a770676bcd3eead3d5fa771c74ee03aca

you can find a tentative workaround for the "deleting holes" issue, at least for MULTIPOLYGON inputs. It is based on "splitting" the MULTIPOLYGON to polygon editing the multipolygon and then putting it back together exploiting the fact that when "casting" from MULTI to POLY the information about the "origins" of the different single polys is stored in the row.names of the casted object.

It's a bit of a ugly workaround, so I did not do a PR yet, but maybe it could be a starting point (tested on a couple of datasets and seems to work).