polygon type detection - clockwise/anticlockwise

[richardbeare]

Hi, I'm chasing some recommendations. Over at r-spatial/osmplotr we've been working on reliable ways of generating polygons that represent sea and land. The starting point is a set of what osm calls "ways" that represent sections of coastline. In osm the land is on the left of such ways. An osm query includes a bounding box and one would receive all ways that cross or are included in the bounding box. However for plotting purposes we need polygons, if we want to colour the sea/land differently. Our original procedure used non-sf procedure attempted to create polygons just outside the bounding box, and it worked some of the time, but broke if there were too many way/bbox crossings.

using sf functions of line_merge, polygonize etc, I can create a section of land and sea polygons, but I need to tell which is which.

So, my question is:

1. Will anything about the original osm ordering of points get preserved? If it is preserved I ought to be able to use clockwise/anticlockise to separate sea/land
2. Is there a clockwise/anticlockwise classifier for polygons?

Any better approaches?

Thanks

[edzer]

Hi, I guess you meant ropensci/osmplotr? Re: your questions:

1. this is entirely up to the osm reader, you didn't tell whether you used osmdata for this or the OSM driver through sf::st_read; the former one is in your hands, so to speak, the latter is highly configurable (see the driver docs and the osmconf.ini file)
2. not a classifier, but since sf 0.5-6 there is a (non-exported) check_ring_dir function that corrects wrongly directed rings; it is also an argument to st_read, st_sf and st_sfc and reverts rings if they have unexpected direction (simple features expect counter-clockwise outer rings, clockwise inner rings)

I guess your real problem is that you don't have rings (polygons) in the first place.

[rsbivand]

Any relationship to this on the GEOS list?

[richardbeare]

The data is coming in via a query from the R osmdata package.

I'm not aware of any direct connection to that GEOS list post.

The problem is that coastline needs to be assembled into connected lines, then combined with the binding box. The following trial code gets most of the way there, but there are two polygons for each, so I need to check which is which. I'll look into the check_ring_direction to see if it helps.

sf_sea <- function(osmd, bbox) {
  ## osmd is a coastline query, and could contain
  ## both lines and polygons.
  ## Lines need to be merged and converted to polygons.
  ## Both need to be clipped against the bounding box
  bbxcoords <- rbind(c(bbox[1, 2], bbox[2, 2]),
                     c(bbox[1, 2], bbox[2, 1]),
                     c(bbox[1, 1], bbox[2, 1]),
                     c(bbox[1, 1], bbox[2, 2])
  )
  bbxcoords <- rbind(bbxcoords, bbxcoords[1,])
  bbx <- st_sfc(st_cast(st_multipoint(bbxcoords), "POLYGON"), crs=st_crs(osmd$osm_lines))
  ## merge lines
  if (!is.null(osmd$osm_lines)) {
    m1 <- st_cast(st_line_merge(st_union(osmd$osm_lines)), "LINESTRING")
    m2 <- st_polygonize(m1)
    k <- st_dimension(m2)
    ## Cast to polygon - need to check what happens when there aren't any
    m2 <- st_cast(m2)
    ## collect the ones remaining as lines
    m1 <- m1[is.na(k)]
    ## clips off the outside parts of line
    m1 <- st_intersection(bbx, m1)
    ## union with the lines in the polygon
    bbxl <- st_cast(bbx, "LINESTRING")
    m1 <- st_union(m1, bbxl, by_feature=TRUE)
    polys <- st_cast(st_polygonize(m1))
    ## now to figure out which polygons are sea
    ## sf has outer 
    ## Original osmdata will have land on the left.

  }
  ## intersect box and lines
  return(list(bbx, polys))
}

[edzer]

@rsbivand does GEOSNormalize_r correct ring directions? I looked for this, but couldn't find it. Are there other docs to the c api then the geos_c.h header file?

@richardbeare you may want to read your way into what a topological representation of a choropleth map means, and how you go from there to ring representations. Possibly, the OSM GDAL driver could help with this. Examples with code would be more helpful then code snippets, but in the end, this seems a problem of osmdata, since simple features don't care about topology.

[richardbeare]

Apologies - I might not have been clear. I don't think there are problems in how sf handles topology - more after advice on a particular problem associated with osm data, and I think that the way osm structures this data is an essential part of solving it. Here's some code, for the record:

library (osmplotr)
library (osmdata)
library (magrittr)

bbox <- osmdata::getbb ("greater melbourne, australia")
coast <- opq (bbox = bbox) %>%
    add_osm_feature (key = "natural", value = "coastline") %>%
    osmdata_sf (quiet = FALSE)

## using github version of ggplot2:
ggplot() + geom_sf(data=coast$osm_lines)

Note that there are also some polygons, that we haven't plotted, and that the island is described by lines ("ways" in osm). Thus the island isn't currently represented as a polygon. Also note that the coastline is broken because the bounding box cuts through a peninsula.

I'd like to be able to fill the land and/or sea, which requires generating a polygon form.

Using sf_sea, above

pp <- sf_sea(coast, bbox)
pp[[2]]

Geometry set for 4 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 144.444 ymin: -38.49937 xmax: 146.1925 ymax: -37.40175
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POLYGON ((144.44405 -38.2824866548647, 144.4444...
POLYGON ((144.44405 -38.150186530243, 144.44420...
POLYGON ((144.44405 -38.0913574222288, 144.4442...
POLYGON ((145.429836636393 -38.49937, 145.42986...

Now, if we plot the first two polygons, we see that they contain the peninsula, with different parts of the bounding box completing the polygon. If there was a way of telling which was sea and which was land, we'd be well on the way to solving the bigger problem

plot(pp[[2]][1])
plot(pp[[2]][2])

Currently having trouble with gdal on the mac, so will need to test under linux soon.

Hope this makes the problem clearer - definitely not suggesting sf bugs.

[edzer]

Thanks, much clearer, and it makes it also clearer that this is not an issue of sf. If your edges tell you what is left and what is right, use that. Maybe post on SO?

[edzer]

@rsbivand normalize described here - rings clockwise, holes counter clockwise - whatever ...