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

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.


 Jeff

	[[alternative HTML version deleted]]

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.


 Jeff

	[[alternative HTML version deleted]]

I suggest taking this question to r-sig-geo, if you haven't already.
-Don

-- 
Don MacQueen

Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062





On 9/18/13 9:07 PM, "Jeff Marcus" <jeff.n.marcus at gmail.com>
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.
>
>
> Jeff
>
>	[[alternative HTML version deleted]]
>
>______________________________________________
>R-help at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide
>http://www.R-project.org/posting-guide.html
>and provide commented, minimal, self-contained, reproducible code.

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.