(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.biv...@nhh.no
_______________________________________________
geos-devel mailing list
geos-devel@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/geos-devel
