[R-sig-Geo] removing internal boundaries from a spatstat window that was created from three adjacent counties

Christopher W. Ryan cryan at binghamton.edu
Mon Feb 29 19:43:20 CET 2016

Adrian, Rolf, et al--

Thanks so much; worked like a charm. Any value of r below 1 (meter, I
assume?) seemed to serve well.

Missing FAQ aside, your book is phenomenal. I'm working through it in
concert with my current project and learning a ton.


Christopher W. Ryan, MD, MS

Early success is a terrible teacher. You’re essentially being rewarded
for a lack of preparation, so when you find yourself in a situation
where you must prepare, you can’t do it. You don’t know how.
--Chris Hadfield, An Astronaut's Guide to Life on Earth

Adrian Baddeley wrote:
> Dear Christopher
> Try
>          newW <- closing(sremsWindow, r) 
> where 'r' is a very small distance (slightly greater than the largest gap that you want to rub out).
>  See closing.owin
> (Dang, that should have been a FAQ in our book)
> A
> Prof Adrian Baddeley DSc FAA
> Department of Mathematics and Statistics
> Curtin University, Perth, Western Australia
> ________________________________________
> From: Rolf Turner <r.turner at auckland.ac.nz>
> Sent: Tuesday, 16 February 2016 6:03 AM
> To: Christopher W. Ryan
> Cc: r-sig-geo at r-project.org; Adrian Baddeley; Ege Rubak
> Subject: Re: [R-sig-Geo] removing internal boundaries from a spatstat window that was created from three adjacent counties
> On 13/02/16 11:59, Christopher W. Ryan wrote:
>> In spatstat, I want to create a single window from three adjacent
>> counties in New York State, USA.  My original data is a shapefile,
>> "cty036" showing all the New York State county boundaries.  Here's what
>> I've done so far:
>> ========== begin code =========
>> library(shapefiles)
>> library(spatstat)
>> library(maptools)
>> library(sp)
>> library(RColorBrewer)
>> library(rgdal)
>> nysArea <- readOGR("C:/DATA/GeographicData", "cty036")
>> nysAreaUTM <- spTransform(nysArea, CRS("+proj=utm +zone=18 +ellps=GRS80
>> +units=m +no_defs"))
>> sremsAreaUTM <- subset(nysAreaUTM, NAME=="Broome" | NAME=="Tioga" |
>> NAME=="Chenango")
>> sremsWindow <- as(sremsAreaUTM, "owin")
>> ======== end code =============
>> This works pretty well, producing a 3-county owin object. But the
>> dividing lines between the counties are shown in the window, whenever I
>> plot it or plot subsequent point patterns that use the window.  In
>> essence, in plots it looks like 3 polygons instead of one big one. I'd
>> prefer not to have the inter-county boundaries be visible--I'd rather
>> have just one big polygon for the whole area. How can I remove them?  Or
>> should I create the window differently from the get-go?
>> Thanks.
> I think it probable that your owin object has come out as a multiple
> polygon.  Look at length(sremsWindow$bdry) --- this is probably 3 (or more).
> I would guess that the internal boundaries actually consist of *pairs*
> of closely juxtaposed lines --- which look like single lines when plotted.
> Have you read the spatstat vignette "shapefiles"?  That might give you
> some guidance as to how to proceed.
> I would have thought the automatic repair process that spatstat now uses
> would fix up this problem.  What version of spatstat are you using?
> You *might* be able to fix things up by doing (something like):
> W1 <- owin(poly=sremsWindow$bdry[[1]])
> W2 <- owin(poly=sremsWindow$bdry[[2]])
> W3 <- owin(poly=sremsWindow$bdry[[3]])
> W  <- union.owin(W1,W2,W3)
> But if my guess about the internal boundaries actually being pairs of
> line is correct, this may not work.
> It's hard to say without having your sremsWindow object to experiment on.
> Or perhaps you need to extract the three counties separately as owin
> objects and the apply union.owin() to the three results.  (Again, if my
> guess is correct, this may not work.)
> Adrian or Ege may be able to propose a more efficient/effective solution.
> cheers,
> Rolf Turner
> --
> Technical Editor ANZJS
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276

More information about the R-sig-Geo mailing list