[R-sig-Geo] [R-sig-geo] How do I ensure that the polygon in spatstat::owin(poly=<polygon>) does not have ?negative area?

[Adrian Baddeley]

Jeff Marcus <jeff.n.marcus at gmail.com> writes:

> I am a new user of the R spatstat package and am having problems creating a
> polygonal observation window with owin(). 

   [ ...]

>  Error in if (w.area < 0) stop(paste("Area of polygon is negative -",
> "maybe traversed in >wrong direction?")) : missing value where TRUE/FALSE
> needed

If you look carefully at this report, it does not say that the area of the polygon was negative.
It says:
    Error in <<<some code I was trying to execute>>>:  missing value where TRUE/FALSE needed

The only variable in the <<code>> is 'w.area', so this report means that the value of w.area was NA or NULL.

That's not very informative unless you go drilling into the spatstat code, but I can tell you what it means: 
when spatstat tried to calculate the polygon area to check whether it was positive or negative, 
the calculated area 'w.area' was NA or NULL.

So ... Are there any NA's in your polygon coordinates? 

(Spatstat does not currently check for NA's - I will fix that!)

Can you plot the coordinates using polygon()? 

If everything seems ok, please send me the data and I'll try to figure out what is happening.

regards
Adrian Baddeley - spatstat coauthor

Prof Adrian Baddeley FAA
University of Western Australia  /and/ CSIRO Mathematics, Informatics & Statistics
Mail: <cet.uwa.edu.au>             Skype: adrian.baddeley

------------------------------

Message: 3
Date: Thu, 19 Sep 2013 16:45:52 -0400
From: 
To: r-sig-geo at r-project.org
Subject: [R-sig-Geo] How do I ensure that the polygon in
        spatstat::owin(poly=<polygon>) does not have ?negative area?
Message-ID:
        <CAHAP1TC=jqmhvdRYDKCtCkXv30jk2H5Q51cuksizqi9cj8kG2A at mail.gmail.com>
Content-Type: text/plain

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]]