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