How do I ensure that the polygon in spatstat::owin(poly=) does not have "negative area"

Question

I am a new user of the R spatstat package and am having problems creating a polygonal observation window with owin(). Code follows:

library("maps")
library ("sp")` 
library("spatstat")
mass.map <- map("state", "massachusetts:main", fill=T) # This returns a data frame includding x and y components that form a polygon of massachusetts mainland`

mass.win <- owin(poly=data.frame(x=mass.map$x, y=mass.map$y)

Error in if (w.area < 0) stop(paste("Area of polygon is negative -", "maybe traversed in >wrong direction?")) : missing value where TRUE/FALSE needed

I tried things like reversing the order of the polygon and got same error.

 mass.win <- owin(poly=data.frame(x=rev(mass.map$x), y=rev(mass.map$y)))

Polygon contains duplicated vertices

Polygon is self-intersecting Error in owin(poly = data.frame(x = rev(mass.map$x), y = rev(mass.map$y))) : Polygon data contain duplicated vertices and self-intersection

Then I figured that maybe the polygon returned by map() is not meant to be fed to owin(). So I tried loading a massachusetts shape file (I am totally taking guesses at this point).:

x <- readShapePoly("../Geog/OUTLINE25K_POLY") ## The shape file for MASS, loaded from MassGIS website
mass.poly <- x <- readShapePoly("../Geog/OUTLINE25K_POLY", force_ring=T, delete_null_obj=T) ## I got following error whether or not I used force_ring

mass.owin <- as(mass.poly, "owin") Checking 1006 polygons...1, Polygon 1 contains duplicated vertices [Checking polygon with 91844 edges...] 2, 3, .. [etd 1:21:52] ....10 [etd 36:12] ..... [etd 23:10] ....20 [etd 16:59] ..... [etd 13:22] ....30 [etd 11:01] ..... [etd 9:21] ....40 [etd 8:06] ..... [etd 7:09] ....50 [etd 6:23] ..... [etd 5:46] ....60 [etd 5:15] ...[Checking polygon with 2449 edges...] .. [etd 4:49] ....70 [etd 4:27] ..... [etd 4:07] ....80 [etd 3:50] ..... [etd 3:36] ....90 [etd 3:22] ..... [etd 3:11] ....100 [ etc.

I got messages complaining about intersecting vertices, etc. and it failed to build the polygon.

Some context on problem: I am trying to use functions in spatstat for spatial relative risk calculations, i.e, the spatial ratio of denstity of cases vs. controls. For that I need an observation window and point plot within that window. I could cheat and make the observation window a rectangle around massachusetts but that would presumably distort values near the coast. In any case, I'd like to learn how to do this right for any future work I do with this package. Thanks for any help you can provide.

Answer 1

EDIT 2021-02-04: I had this clockwise/anti-clockwise problem again and found my own answer the only one addressing the issue for spatstat. The answer was not very helpful. Have improved it to make it useful for a wider range of applications.

spatstat package wants owin object coordinates specified anticlockwise. A quote from the documentation of version 1.36-0:

single polygon: If poly is a matrix or data frame with two columns, or a structure with two component vectors x and y of equal length, then these values are interpreted as the cartesian coordinates of the vertices of a polygon circumscribing the window. The vertices must be listed anticlockwise. No vertex should be repeated (i.e. do not repeat the first vertex).

library("maps")
library ("sp")
library("spatstat")
mass.map <- map("state", "massachusetts:main", fill=T)

First, you'll need to find out whether your polygon is running clockwise or anticlock-wise. The approach here can be used to find out the answer. Formulate it as a function:

#' @title Check whether points for an owin are clockwise
#' @param x a dataframe with x coordinates in the first column and y coordinates in the second. 
#' @details Similarly to owin, the polygon should not be closed
#' @return A logical telling whether the polygon is arranged clockwise.
#' @author The idea has been scavenged from https://stackoverflow.com/a/1165943/1082004

clockwise <- function(x) {
    
  x.coords <- c(x[[1]], x[[1]][1])
  y.coords <- c(x[[2]], x[[2]][1])
  
  double.area <- sum(sapply(2:length(x.coords), function(i) {
    (x.coords[i] - x.coords[i-1])*(y.coords[i] + y.coords[i-1])
  }))
  
  double.area > 0
} 

Now, see whether the coordinates in mass.map are clockwise:

clockwise(data.frame(x=mass.map$x, y=mass.map$y))
#> TRUE

They are. Turn the coordinates anti-clockwise by using rev function, which reverses vectors for x and y:

mass.win <- owin(poly=data.frame(x=rev(mass.map$x), y=rev(mass.map$y)))
plot(mass.win)

Answer 2

The original question was posted 7 years ago (2014). The spatstat package has changed vastly since then.

The error messages about duplicated vertices and self-intersecting polygons do not occur anymore, because the code now automatically corrects these flaws in the polygon data.

The error message Area of polygon is negative - maybe traversed in wrong direction? still occurs. This is because owin requires polygon vertices to be listed in anticlockwise order around an outer boundary (and clockwise order around the boundary of a hole). If your window is a single polygon rather than a collection of polygons, the fix is to simply reverse the order of the coordinate vectors using rev:

library(maps)
library(sp)
library(spatstat)
mass.map <- map("state", "massachusetts:main", fill=TRUE)
mass.win <- owin(poly=lapply(mass.map, rev))

We can't make this happen automatically because it is not always what the user wanted.

For more information, see Chapter 3 of the spatstat book, which you can download for free from here.

Answer 3

I'm quite familiar with the problem you have. If I understood right, you desire to get a polygonal object of class owin. Following codes will help you to acquire what you wanted.

#Loading the required packages
library (sp)
library(maps)

# Reading the dataset
mass.map <- map("state", "massachusetts:main", fill=T)

IDs <- sapply(strsplit(mass.map$names, ":"), function(x) x[1])
# converting can be done by loading the powerful package "maptools"
library(maptools)
# Converting into SpatialPolygons-object
mass.poly <- map2SpatialPolygons(mass.map, IDs = IDs)

# Converting by coercing into owin-object
mass.owin <- as.owin.SpatialPolygons(mass.poly)