[R] How do I ensure that the polygon in spatstat::owin(poly=<polygon>) does not have “negative area”

[Rolf Turner]

You were nearly there.  No need to muck about with the shapefile 
business which
introduces a whole lot of sub-polygons (counties or townships I guess).

The polygon for the border of Massachusetts provided in "maps" has its 
vertices
in clockwise order and its first and last vertices identical to each 
other.  You need
to (a) reverse the direction (as you have done) and (b) remove the 
duplication.
Easiest to do (b) first:

mass.df <- with(mass.map,data.frame(x=rev(x),y=rev(y)))
mass.df <- unique(mass.df)
mass.win  <- owin(poly=mass.df)
plot(mass.win) # OMMMMMMMMM!

     cheers,

     Rolf Turner

P.S.  It is inadvisable to use "T" when you mean "TRUE".  Write out 
"TRUE" in
full.  (The name "T" can be over-written, leading to dangers. E.g. "T <- 
FALSE" !!!
The name "TRUE" cannot be over-written.)

     R. T.

On 09/19/13 16:07, Jeff Marcus wrote:
> 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.
>
> Note: I cross-posted this to STack Overflow and then realized that r-help
> is probably a better forum.