[geos-devel] rgeos interface to R classes (sp, others)
Roger Bivand
Roger.Bivand at nhh.no
Tue May 4 05:07:59 EDT 2010
On Wed, 3 Feb 2010, Roger Bivand wrote:
> (reverting to list - link to test C code below).
For closure on the thread of touching polygons not being merged by union
operations:
The rgeos R/GEOS interface is now a GSoC project within the R project, and
the "student", Colin Rundel, has found the root of the problem: when
creating GEOSCoordSeq, the val values had not been java_math_round'ed.
Once they are java_math_round'ed, rather than just taken as the doubles
they were, the problem goes away. Reading from WKT, etc., couldn't
reproduce the problem, because the read values were rounded on their way
in.
So I wasn't using the default precision model, just passing doubles from
R. Thanks to Colin for his analysis and resolution! Is there a suitable
place to let others know that double val being fed say to
GEOSCoordSeq_setX_r in the C API should apparently be pushed through a
precision filter first?
Roger
>
> 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