editing + deleting in editFeatures does't work

[tim-salabim]

When I do

editFeatures(franconia[1:5, ]) %>% 
  mapview()

and in that session both delete a feature and edit another one, I get the error:

Error in CPL_get_bbox(obj, 3) : Not a matrix.

traceback:

22.
stop(structure(list(message = "Not a matrix.", call = CPL_get_bbox(obj, 
    3), cppstack = NULL), .Names = c("message", "call", "cppstack"
), class = c("Rcpp::not_a_matrix", "C++Error", "error", "condition"
))) 
21.
CPL_get_bbox(obj, 3) 
20.
structure(CPL_get_bbox(obj, 3), names = c("xmin", "ymin", "xmax", 
    "ymax")) 
19.
bbox.MtrxSetSetSet(obj) 
18.
structure(bb, class = "bbox", crs = st_crs(x)) 
17.
bb_wrap(bbox.MtrxSetSetSet(obj), obj) 
16.
st_bbox.sfc_MULTIPOLYGON(x) 
15.
st_bbox(x) 
14.
`[.sfc`(geom, i) 
13.
geom[i] 
12.
`[.sf`(orig, which(!(orig_ids %in% del_ids)), ) at merge.R#72
11.
orig[which(!(orig_ids %in% del_ids)), ] at merge.R#72
10.
merge_delete(structure(list(NUTS_ID = c("DE241", "DE242", "DE243", 
"DE244", "DE245"), SHAPE_AREA = c(0.00673601232111, 0.00842446859281, 
0.00598234136934, 0.00732948035553, 0.146698316001), SHAPE_LEN = c(0.392622499474, 
0.624726298357, 0.518547080489, 0.456981511253, 3.4819699025),  ... 
9.
eval(call(paste0("merge_", op), left_sf, sf_merge, c(edit_id = "layerId"))) 
8.
eval(call(paste0("merge_", op), left_sf, sf_merge, c(edit_id = "layerId"))) at edit.R#199
7.
f(init, x[[i]]) 
6.
Reduce(function(left_sf, op) {
    op <- tolower(op)
    if (op == "add") 
        sf_merge <- crud$finished ... at edit.R#172
5.
editFeatures.sf(franconia[1:5, ]) at edit.R#119
4.
editFeatures(franconia[1:5, ]) 
3.
eval(lhs, parent, parent) 
2.
eval(lhs, parent, parent) 
1.
editFeatures(franconia[1:5, ]) %>% mapview() 

[edzer]

That comes from sf; I didn't debug, but a reason might be that a POLYGON object has class MULTIPOLYGON, or vice versa, or something similar.

[timelyportfolio]

Thanks @tim-salabim for the report, and thanks @edzer for confirming my suspicion. I think this relates to #32.

[edzer]

Doing sth similar, I see

Error in MtrxSetSet(x, dim, type = "MULTIPOLYGON", needClosed = TRUE) : 
  polygons not (all) closed

details:

Details

 traceback()
32: stop("polygons not (all) closed")
31: MtrxSetSet(x, dim, type = "MULTIPOLYGON", needClosed = TRUE)
30: st_multipolygon(list(unclass(x)), dim = dim)
29: POLYGON2MULTIPOLYGON(x)
28: FUN(X[[i]], ...)
27: lapply(x, function(x) if (inherits(x, "POLYGON")) POLYGON2MULTIPOLYGON(x) else x)
26: st_cast_sfc_default(x)
25: st_cast.sfc(st_geometry(x), to, group_or_split = do_split)
24: st_cast(st_geometry(x), to, group_or_split = do_split)
23: st_cast.sf(x)
22: sf::st_cast(x)
21: sf::st_geometry(sf::st_cast(x))
20: getGeometryType(x)
19: standardColor(x)
18: vectorColors(x = x, zcol = zcol, colors = color, at = at, na.color = na.color)
17: leaflet_sf(x, map = map, zcol = zcol, color = color, col.regions = col.regions, 
        at = at, na.color = na.color, cex = cex, lwd = lwd, alpha = alpha, 
        alpha.regions = alpha.regions, na.alpha = na.alpha, map.types = map.types, 
        verbose = verbose, popup = popup, layer.name = layer.name, 
        label = label, legend = legend, legend.opacity = legend.opacity, 
        homebutton = homebutton, native.crs = native.crs, highlight = highlight, 
        maxpoints = maxpoints, ...)
16: .local(x, ...)
15: mapView(...)
14: mapView(...)
13: .Method(...)
12: eval(.dotsCall, env)
11: eval(.dotsCall, env)
10: standardGeneric("mapview")
9: mapview(.)
8: function_list[[k]](value)
7: withVisible(function_list[[k]](value))
6: freduce(value, `_function_list`)
5: `_fseq`(`_lhs`)
4: eval(quote(`_fseq`(`_lhs`)), env, env)
3: eval(quote(`_fseq`(`_lhs`)), env, env)
2: withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
1: editFeatures(franconia[1:5, ]) %>% mapview()
> 

[tim-salabim]

@edzer what did you do? I cannot seem to reproduce the error you see.

[edzer]

Are you using the GH version of sf?

[tim-salabim]

Yes, still the only error I get is Error in CPL_get_bbox(obj, 3) : Not a matrix.

[tim-salabim]

I can now replicate the polygons not (all) closed error. I think it only happens when you edit MULTIPOLYGONs (i.e. the big one or the one east of it. Maybe we need to cast them to POLYGON first, do the editing and then cast them back to MULTIPOLYGON before returning the object? Obviously, this only works for similar type collections (GEOMETRYs).

[tim-salabim]

A quick trial seems to back my suspicion.

editFeatures(st_cast(franconia[1:5, ], "POLYGON")) %>%
  mapview()

does not seem to produce the polygons not (all) closed error.

[timelyportfolio]

Darn, that is a troubling regression. I will research thoroughly tomorrow and hopefully at least trace the source.

[timelyportfolio]

@edzer, @tim-salabim ok, I am going to need some help here to insure the correct behavior. After quite a bit of research, I don't think this is a regression in mapedit. Leaflet.Draw does not handle MULTIPOLYGON, so when we edit franconia, the edited features are returned as POLYGON GeoJSON. So, to reproduce the error without mapedit, we can simply

fran2 <- franconia
st_geometry(fran2[1,]) <- st_geometry(st_cast(franconia[1,],"POLYGON"))
plot(fran2)

which gives us

Error in CPL_get_bbox(obj, 3) : Not a matrix.

Compounding the Problem

In the simple case above with length==1, w could solve by attempting to cast or convert back to the original type before merging. Unfortunately, this does not work if we consider franconia[2,]
with length==2, which Leaflet.Draw returns as one feature of type POLYGON with two sets of coordinates. If I understand correctly, this will only happen with holes, or am I wrong?

If you tell me, this is ok, then we can just do as above and st_cast before merging. If it is not ok, it seems like we might have to make the decision to support holes, which we know Leaflet.Draw does not handle, or ignore holes and convert each set of coordinates to a polygon and make it a MULTIPOLYGON.

[edzer]

As we found out earlier, the

st_geometry(fran2[1,]) <- st_geometry(st_cast(franconia[1,],"POLYGON"))

causes trouble because it has a nested replacement function ([<- and st_geometry<-)

Why don't you simply use

st_cast(franconia, "POLYGON")

to cast franconia to POLYGON? OK, it changes franconia from 37 to 40 features, but isn't that what you want? In case you'd want to merge some of them together later on,

> lengths(st_geometry(franconia))
 [1] 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

points you to the multi ones.

[timelyportfolio]

That is where I need help from more experienced geo people in trying to determine the best "solution". I unfortunately don't have the domain expertise to know if cast to polygon first is acceptable. @tim-salabim, do you have a recommendation?

@edzer, I just can't stop myself with the nested assignment, but that does not seem to be the problem here. Any alternative gives us the same error, since it seems the problem is POLYGON amidst MULTIPOLYGON in sf.

library(mapview)
library(mapedit)
library(sf)

fran2 <- franconia
st_geometry(fran2)[1] <- st_geometry(st_cast(franconia[1,],"POLYGON"))
plot(fran2)

fran3 <- franconia
st_geometry(fran3[1,]) <- st_geometry(st_cast(franconia[1,],"POLYGON"))
plot(fran3)

fran_g <- st_geometry(franconia)
fran_g[1] <- st_geometry(st_cast(franconia[1,],"POLYGON"))
plot(fran_g)

[timelyportfolio]

@edzer, would the correct way to do this be this or a derivative?

fran_g <- st_geometry(franconia)
fran_g[1] <- st_geometry(st_cast(franconia[1,],"POLYGON"))
fran_g2 <- st_sfc(fran_g)
plot(fran_g2)