[R-sig-Geo] distmap for polygons imported from shapefiles

[Agustin Lobo]

The distmap calculation works fine for polygons
imported from shapefiles, thanks Adrian and Roger!
Agus

Roger Bivand escribió:
> 
> Doing all the conversions as specified, it does work:
> 
> library(maptools)
> library(spatstat)
> load("refG.rda") 

#Note, refG comes from:
#refG <- 
#readOGR("C:/ALOBO/dipu2006/espaisoberts_P037/GARRAF/GFO_EOA_TOTS_v4",layer="GFO_EOA_Acti_v4")

> class(refG)
> ?as.owin.SpatialPolygons
> 
> # (note no method for SpatialPolygonsDataFrame)
> 
> plot(refG)
> .spatstat_check <- FALSE # this is the trick used to pass FALSE to check=
> w2 <- as(as(refG, "SpatialPolygons"), "owin")
> rm(.spatstat_check)
> plot(w2)
> p1 <- as.psp(w2)
> plot(p1)
> d <- distmap(p1)
> summary(d)
> plot(d)
> d[w2] <- NA
> summary(d)
> plot(d)
> 
> Use the ... to distmap() to set your own raster parameters.
> 
> Hope this helps,
> 
> Roger
> 
> 
>> I get:
>>
>> Warning messages:
>> 1: $ operator is deprecated for atomic vectors, returning NULL in: po$x
>> 2: $ operator is deprecated for atomic vectors, returning NULL in: po$x
>> 3: $ operator is deprecated for atomic vectors, returning NULL in: po$y
>> 4: $ operator is deprecated for atomic vectors, returning NULL in: po$x
>> 5: $ operator is deprecated for atomic vectors, returning NULL in: po$y
>>
>> If I represent
>>
>>>  plot(w)
>> I get an empty map. (summary(w) does not work)
>> and also for d:
>>>  summary(d)
>> real-valued pixel image
>> 100 x 100 pixel array (ny, nx)
>> enclosing rectangle: [8.67361737988404e-19, 1] x 
>> [8.67361737988404e-19, 1] units
>> dimensions of each pixel: 0.01 x 0.01 units
>> Image is defined on the full rectangular grid
>> Frame area =  1 square units
>> Pixel values :
>>         range = [0,0]
>>         integral = 0
>>         mean = 0
>>
>> While
>>>  plot(refG)
>> is correct.
>>
>> Finally,
>>>  d[w] <- NA
>> Error in verify.xypolygon(polly) : polygon must be a list with two 
>> components x and y
>>
>> I attach refG again. Remember it comes from
>> refG <- 
>> readOGR("C:/ALOBO/dipu2006/espaisoberts_P037/GARRAF/GFO_EOA_TOTS_v4",layer="GFO_EOA_Acti_v4") 
>>
>>
>> Agus
>>
>>
>>
>> adrian at maths.uwa.edu.au escribió:
>>>  Agustin and Roger -
>>>
>>>  Regarding your query about the distance map of polygonal windows:
>>>
>>>  I forgot to mention that you can compute the distance map of a 
>>> collection
>>>  of *line segments* using distmap.psp. This is computed exactly using
>>>  Euclidean geometry. You can convert a polygonal window to a psp object
>>>  using as.psp.
>>>
>>>  So the following would work safely. Suppose 'mypoly' is a collection of
>>>  polygons that do not pass the polygon checking code in spatstat because
>>>  they have some common vertices. Then
>>>                    w <- owin(poly=mypoly, check=FALSE)
>>>                    d <- distmap(as.psp(w))
>>>                    d[w] <- NA
>>>
>>>  I hope this helps
>>>
>>>  regards
>>>  Adrian
>>>
>>>
>>>
>>
>>
> 

-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster