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@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.
