[geos-devel] rgeos interface to R classes (sp, others)
Roger Bivand
Roger.Bivand at nhh.no
Wed Feb 3 13:28:50 EST 2010
(reverting to list - link to test C code below).
On Tue, 2 Feb 2010, Martin Davis wrote:
> Roger,
>>
>>
>> Thanks for checking. In my coding using GEOS C API they do not, leading to
>> my puzzlement. I guess it isn't GEOS, so must be my coding, but I can't see
>> where - nothing I do (I think) is scale-dependent.
> Very strange. It sounds like it might be a Precision Model problem - but I'm
> guessing you're just using the default precision model, which is floating
> point and should handle this test case fine.
>
> GEOS gets extremely well tested via its heritage in JTS and by its use in
> PostGIS, so the problem is almost certainly in how it is being called, or in
> some issue with the input data. I suggest printing out the input and results
> immediately before and after the GEOS calls, to see if something is going
> wrong there.
>
> I've tried looking at your code a bit more, but I have to confess that I get
> lost in the mixture of function calls to both R and GEOS. Perhaps you could
> try and make a simple test program using only GEOS calls? That will be
> needed in any case if there is truly a GEOS bug.
Martin:
I've put a small test program on:
http://spatial.nhh.no/misc/test.c
which doesn't reproduce the problem exactly, but obviously isn't right. It
may be a validity problem, with:
Duplicate Rings at or near point 0 0
type MultiPolygon, valid: 0
reason: Duplicate Rings[0 0]
in (now) both cases, even though I don't see that the three polygons in
the multipolygon duplicate anything. I now get Polygon objects that are
not disjoint and have 0 distance not touching. Here UnionCascaded also
breaks down (probably because of the validity problem). The data are
hard-coded (and may be wrong too?), rather than brought in from WKT.
Roger
>>
>>>
>>> But perhaps your problem lies elsewhere? I don't understand your
>>> reference to GEOMTouches in the first email - perhaps this is not
>>> producing the right result, or is being used incorrectly?
>>
>> Trying to find out if it was my coding, I tried sequential GEOMTouches
>> among the polygons in the multipolygon. Again, with a 0.1 x 0.1 square,
>> they do not touch, but with 100 x 100 they do. I found the problem after
>> giving advice on the R-sig-geo list about how to create a vector patch map
>> from a raster (dissolve borders between contiguous vectorised cells), and
>> was most surprised to see the toy example fail.
> Ah, ok, I understand.
>>
>> The underlying motivation is to provide R/GEOS integration without Java -
>> one contributed R package (cshapes) does bundle JTS through Java, but
>> exchanges data by writing and reading shapefiles (R doesn't fit Java well
>> installation-wise). In addition, the port of GPC to R that I've used for
>> this before has serious license problems, and JTS/GEOS are both more
>> general and license-wise clean.
> That makes perfect sense. It's always nicer to have a library that is native
> to the platform you're working on. (I'm amazed that someone actually
> interfaced to JTS via files and presumably system calls - but I guess when
> all you have is a hammer...)
>
>
>
>>
>> Roger
>>
>>>
>>>
>>> Roger Bivand wrote:
>>>> On Mon, 1 Feb 2010, Martin Davis wrote:
>>>>
>>>>> Can you post the WKT/WKB of your test case?
>>>>
>>>> If this is WKT, then:
>>>>
>>>>> SP2WKT(spol0)
>>>> [1] "MULTIPOLYGON(((0 0.1,0 0.2,0.1 0.2,0.1 0.1,0 0.1)),
>>>> ((0.1 0.1,0.1 0.2,0.2 0.2,0.2 0.1,0.1 0.1)),
>>>> ((0 0,0 0.1,0.1 0.1,0.1 0,0 0)))"
>>>>> SP2WKT(spolx0)
>>>> [1] "MULTIPOLYGON(((0 100,0 200,100 200,100 100,0 100)),
>>>> ((100 100,100 200,200 200,200 100,100 100)),
>>>> ((0 0,0 100,100 100,100 0,0 0)))"
>>>>>
>>>>
>>>> I don't think that there is an OGR driver, is there? So this is a text
>>>> representation for the two objects, each with 3 Polygon members, and only
>>>> differing in scale (the coordinates of the second are 1000 times those of
>>>> the first). I see the same difference on GEOS 3.2.0/C API 1.6.0 on x86_64
>>>> RHEL5 and with an MSYS-built GEOS 3.2.0/C API 1.6.0 on Win 32 R.
>>>>
>>>> Hope this is good enough, grateful for any ideas,
>>>>
>>>> Roger
>>>>
>>>>
>>>>> Roger Bivand wrote:
>>>>>> On Mon, 1 Feb 2010, Martin Davis wrote:
>>>>>>
>>>>>>> Roger,
>>>>>>>
>>>>>>> The technique of using buffer(0) to union polygons is now deprecated
>>>>>>> in favour of using Unary Union (Geometry.union() - not sure what the
>>>>>>> exact GEOS signature is). Unary Union is usually faster and more
>>>>>>> robust than the previous technique. You might want to check this out.
>>>>>>
>>>>>> Martin,
>>>>>>
>>>>>> Thanks, I've replaced buffer(0) by GEOSUnionCascaded (C API). But the
>>>>>> artefact remains. For input MULTIPOLYGON:
>>>>>>
>>>>>> ---------
>>>>>> | | |
>>>>>> | | |
>>>>>> |--------
>>>>>> | |
>>>>>> | |
>>>>>> -----
>>>>>>
>>>>>> I still get the same MultiPolygon out when the bounding box is
>>>>>> (0,0)(0.2,0.2), but correctly Polygon:
>>>>>>
>>>>>> ---------
>>>>>> | |
>>>>>> | |
>>>>>> | -----
>>>>>> | |
>>>>>> | |
>>>>>> -----
>>>>>>
>>>>>> when the bounding box is (0,0)(200,200). May this be related to the
>>>>>> input Polygon objects only having 5 coordinates (4 corners and
>>>>>> closure)? The generating code is inserting exactly the same doubles
>>>>>> into the coordinates.
>>>>>>
>>>>>> Anyway, I'm sure that Unary Union is a better solution than buffer(0)
>>>>>> in general, so thanks for that!
>>>>>>
>>>>>> Roger
>>>>>>
>>>>>>
>>>>>>>
>>>>>>> Martin
>>>>>>>
>>>>>>> Roger Bivand wrote:
>>>>>>>> Hi,
>>>>>>>>
>>>>>>>> I've felt that I've been making reasonable progress with interfacing
>>>>>>>> GEOS geometries and methods for the R language (cran.r-project.org),
>>>>>>>> in a draft contributed package rgeos:
>>>>>>>>
>>>>>>>> https://r-forge.r-project.org/projects/rgeos/
>>>>>>>>
>>>>>>>> However, I encountered a problem that I do not understand, and would
>>>>>>>> be very grateful if I could be put back on track. The specific
>>>>>>>> problem is that an R function uses the C API to dissolve polygon
>>>>>>>> borders for adjacent polygons sharing a category given as its second
>>>>>>>> argument.
>>>>>>>>
>>>>>>>> I've used the buffer technique from the JTS documentation, and all
>>>>>>>> was well until I tried to dissolve borders between touching squares
>>>>>>>> when the coordinate measures were small (square side 0.1). When the
>>>>>>>> squares are 100 units, all is well, and GEOMTouches is TRUE. But for
>>>>>>>> 0.1, GEOMTouches is FALSE, and no dissolving takes place.
>>>>>>>>
>>>>>>>> This can be reproduced (I'm on RHEL5, x86_64) by installing R,
>>>>>>>> contributed packages sp and maptools from CRAN, GEOS (3.1.1 or
>>>>>>>> 3.2.0), and installing rgeos from R-Forge. Start R, say
>>>>>>>>
>>>>>>>> library(rgeos) example(unionSpatialPolygonsGEOS)
>>>>>>>>
>>>>>>>> and look for undissolved borders in:
>>>>>>>>
>>>>>>>> image(grd, axes=TRUE)
>>>>>>>> plot(spol1, add=TRUE)
>>>>>>>>
>>>>>>>> but dissolved in:
>>>>>>>>
>>>>>>>> image(grdx, axes=TRUE)
>>>>>>>> plot(spol1x, add=TRUE)
>>>>>>>>
>>>>>>>> I started on an alternative implementation using GEOMTouches in
>>>>>>>> unionSpatialPolygonsGEOS(..., buffer=FALSE) output to console, where
>>>>>>>> one sees in the example output:
>>>>>>>>
>>>>>>>>> spol1F <- unionSpatialPolygonsGEOS(as(spol, "SpatialPolygons"),
>>>>>>>> + as.character(spol$xvs), buffer=FALSE)
>>>>>>>> # 4 squares within (0,0), (0.2,0.2) NE, NW, SW share category, SE
>>>>>>>> doesn't
>>>>>>>> npls: 3, nnpls: 3
>>>>>>>> type[0] Polygon
>>>>>>>> i: 0, j: 1, touches: 0
>>>>>>>> i: 0, j: 2, touches: 0
>>>>>>>> type[1] Polygon
>>>>>>>> i: 1, j: 2, touches: 0
>>>>>>>> type[2] Polygon
>>>>>>>> out of function
>>>>>>>> npls: 1, nnpls: 1
>>>>>>>> type[0] Polygon
>>>>>>>> out of function
>>>>>>>>> spol1xF <- unionSpatialPolygonsGEOS(as(spolx, "SpatialPolygons"),
>>>>>>>> + as.character(spolx$xvs), buffer=FALSE)
>>>>>>>> # 4 squares within (0,0), (200,200), same categories
>>>>>>>> npls: 3, nnpls: 3
>>>>>>>> type[0] Polygon
>>>>>>>> i: 0, j: 1, touches: 1
>>>>>>>> i: 0, j: 2, touches: 1
>>>>>>>> type[1] Polygon
>>>>>>>> i: 1, j: 2, touches: 1
>>>>>>>> type[2] Polygon
>>>>>>>> out of function
>>>>>>>> npls: 1, nnpls: 1
>>>>>>>> type[0] Polygon
>>>>>>>> out of function
>>>>>>>>
>>>>>>>> The C code is in:
>>>>>>>>
>>>>>>>> https://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/pkg/src/rgeos_sp.c?rev=44&root=rgeos&view=log
>>>>>>>> around line 466 rgeos_SpatialPolygonsUnion(), calling
>>>>>>>> rgeos_plspairUnion() - desperate test framework, or
>>>>>>>> rgeos_plsbufUnion() which I had thought worked, but which clearly
>>>>>>>> doesn't dissolve small squares.
>>>>>>>>
>>>>>>>> Very grateful for any help, this is a show-stopper, and I had hoped
>>>>>>>> to release rgeos before February (small chance now!)
>>>>>>>>
>>>>>>>> Roger
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>
>
>
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no
More information about the geos-devel
mailing list