Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

as.polygons results in clockwise polygons #1249

Closed
srfall opened this issue Aug 1, 2023 · 3 comments
Closed

as.polygons results in clockwise polygons #1249

srfall opened this issue Aug 1, 2023 · 3 comments

Comments

@srfall
Copy link
Contributor

srfall commented Aug 1, 2023

I'm not an expert on the wkt standard but as far as I am aware of, the exterior boundary of a polygon should be defined counter clockwise and many other applications require this format (e.g see the issue in rgbif here: ropensci/rgbif/issues/662).

Observed problem:

terra::geom(terra::as.polygons(terra::ext(2, 18, 45, 55), crs = "+proj=longlat"), wkt = TRUE)

returns
"POLYGON ((2 45, 2 55, 18 55, 18 45, 2 45))"
which is clockwise, instead of the desired counter clockwise
"POLYGON ((2 45, 18 45,18 55, 2 55, 2 45))"

It looks like as.polygons is the cause:

wkt_polygon_input <- "POLYGON ((2 45, 18 45,18 55, 2 55,  2 45))"

terra::geom(terra::as.polygons(terra::vect(wkt_polygon_input, crs = "+proj=longlat")), wkt = TRUE)
# "POLYGON ((2 45, 2 55, 18 55, 18 45, 2 45))"
# undesired result

terra::geom(terra::vect(wkt_polygon_input, crs = "+proj=longlat"), wkt = TRUE)
# "POLYGON ((2 45, 18 45, 18 55, 2 55, 2 45))"
# desired result
@srfall
Copy link
Contributor Author

srfall commented Aug 1, 2023

If this is not a bug, but rather intended, it might be nice to have a parameter that can specify the desired winding order.

@edzer
Copy link
Contributor

edzer commented Aug 5, 2023

Some background on polygon orientation: r-spatial/sf#2096

@rhijmans
Copy link
Member

I have added a new method forceCCW that returns counter-clockwise polygons.

library(terra)
#terra 1.7.42
p <- terra::as.polygons(terra::ext(2, 18, 45, 55), crs = "+proj=longlat")
x <- forceCCW(p)

geom(p, wkt = TRUE)
#[1] "POLYGON ((2 45, 2 55, 18 55, 18 45, 2 45))"
geom(x, wkt = TRUE)
#[1] "POLYGON ((2 45, 18 45, 18 55, 2 55, 2 45))"

#Also see
plot(x, axes=F, border="azure", lwd=4)
text(as.points(x), 1:4)

It would be trivial to change the order returned by as.polygons<SpatExtent> to CCW if that is helpful.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants