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