Polygon hole is filled when the order of vertices is the same as in main polygon

[tsamsonov]

It seems that the order of vertices in a polygon hole influences the way it is plotted. The problem arises when the POLYGON sfg is constructed manually:

main <- matrix(c( # The order is clockwise
  1, 0,
  0, 2,
  2, 3,
  4, 2,
  3, 0.5,
  1, 0
), ncol = 2, byrow = TRUE)

hole <- matrix(c( # The order is clockwise
  2, 1,
  2, 2,
  3, 2,
  3, 1.5,
  2, 1
), ncol = 2, byrow = TRUE)

pol <- st_polygon(list(main, hole))

plot(pol, col = 'grey', main = 'Clockwise-Clockwise') # the hole is filled

st_is_valid(pol) # reports TRUE

hole2 <- matrix(c( # The order is counterclockwise
  2, 1,
  3, 1.5,
  3, 2,
  2, 2,
  2, 1
), ncol = 2, byrow = TRUE)

pol2 <- st_polygon(list(main, hole2))

plot(pol2, col = 'grey', main = 'Clockwise-Counterclockwise') # the hole is not filled

st_is_valid(pol2) # reports TRUE

In both cases st_is_valid() reports that polygon is valid, but in the first case it is not.

I suggest that the rules of point ordering should be documented explicitly in sf package (in correspondence with OGC Simple Feature Access standard). Some behaviour can be programmed when the order of vertices in a hole is the same as in the main polygon: it can be changed automatically or an error is raised by st_polygon(). Alternatively st_is_valid() may give FALSE in case of inconsistent vertex ordering in polygons and their components.

[rsbivand]

You'd need to look at the polypath documentation, as:

par(mfrow=c(1,2))
plot(pol, col = 'grey', main = 'Clockwise-Clockwise, winding', rule="winding")
plot(pol, col = 'grey', main = 'Clockwise-Clockwise, evenodd', rule="evenodd")

does what you want. Conversion to sp works:

library(sp)
sp_pol <- as(st_geometry(pol), "Spatial")
plot(sp_pol, col = 'grey', main = 'Clockwise-Clockwise, winding', rule="winding")
plot(sp_pol, col = 'grey', main = 'Clockwise-Clockwise, evenodd', rule="evenodd")

plot(st_as_sfc(sp_pol), col = 'grey', main = 'Clockwise-Clockwise, winding', rule="winding")
plot(st_as_sfc(sp_pol), col = 'grey', main = 'Clockwise-Clockwise, evenodd', rule="evenodd")
par(mfrow=c(1,1))

So, yes, sp is more cautious about ring direction, and if it knows that a hole is really a hole, will change the direction if inappropriate.

[edzer]

Where in the SFA 1.2 standard document did you read that the validity of a polygon constrains direction of rings? In any case, we simply pass this problem on to the GEOS function GEOSisValidReason_r, which I presume follows the list of issues SFA mentions in the context of what is meant by valid, here.

Recent https://github.com/r-spatial/sf/commit/f9296739899825e6ff36f0b28e65a2d27cca764a adds ring direction checking, and correction, to sf; you'll get the right ring dirs with

st_sfc(pol, check_ring_dir = TRUE)

I didn't implement check_ring_dir at the st_polygon or st_multipolygon level, but we still could (these functions are IMO mostly used for example construction).

I'm not sure we should do this by default; there is an overhead. @Robinlovelace ?

[tsamsonov]

@rsbivand @edzer Thanks for your replies. If the formal definition of validity does not constrain the direction of rings, then st_is_valid() should not be modified. If check_ring_dir() gives a significant overhead then at least the brief guidelines for manual generation of polygons with holes using st_polygon() should be added to sf documentation.

[edzer]

It is, in ?st_read:

check_ring_dir: logical; if TRUE, polygon ring directions are checked
          and if necessary corrected (when seen from above: exterior
          ring counter clockwise, holes clockwise)

[rsbivand]

Maybe even change the default rule and leave check_ring_dir FALSE?

> benchmark(winding=plot(pol, col = 'grey', rule="winding"), evenodd=plot(pol, col = 'grey', rule="evenodd"), replications=10000)
     test replications elapsed relative user.self sys.self user.child sys.child
2 evenodd        10000  25.500    1.000     17.11    0.761          0         0
1 winding        10000  25.666    1.007     17.14    0.781          0         0

[tsamsonov]

@edzer What if the user is not using st_read? Manual generation, for example. May be a short note on ring direction should be added to sfg: simple feature geometry section of the first vignette. Or fifth vignette on plotting simple features.

[rsbivand]

Aren't the effects seen here rather about the R base plot system for polygons?

[tsamsonov]

@rsbivand If the hole ring direction only affects the plotting result and does not manifest itself in other scenarios of using simple features, then you are right. But the idea remains the same: instruct the user how to prepare sf geometry for correct plotting.

[edzer]

I don't think there is such a thing as correct plotting in this case; in any case both base and grid (ggplot2) are now robust against ring misspecs: https://github.com/r-spatial/sf/commit/5bae05d5050986199f9cf8e0d56d470f0479e5ae and https://github.com/r-spatial/sf/commit/889b8c00b4dfd009447f80d0bf904e169af7a5c1

[tsamsonov]

@edzer Thank you.

[rsbivand]

Looks robust now; agree with @edzer that correct plotting isn't clearly defined.

[tsamsonov]

If the hole of a polygon is filled -- then in most of the cases the image is incorrect from geographic/cartographic/perceptional perspective, because the area of a hole does not belong to the area of a polygon. But I agree that this may not make sense from the computational/graphical point of view. And may be there are some rare cases when the hole should be filled (but I cannot come up with an example). Anyway the graphical issue seems to be worked out now.

[edzer]

Thanks both, for the help. I guess this is one of these cases were the SFA standard would have been more helpful, if it were more strict.

[edzer]

Of course, the new winding rule ("evenodd") breaks on overlapping holes, as in

library(sf)
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
plot(st_polygon(list(m * 3.5, m + 1, m + 1.5)), col = 'grey')

but this kind can at least be found by st_is_valid:

> st_is_valid(st_polygon(list(m * 10, m + 1, m + 1.5)))
[1] FALSE
Warning message:
In evalq((function (..., call. = TRUE, immediate. = FALSE, noBreaks. = FALSE,  :
  Self-intersection at or near point 2 1.5