I would be happy to write a polygon clipping algorithm. I am not much of a
GIS expert, though, so I would need some hand holding - particularly for
the design requirements and definitions.
* what functions are needed?
* what should their inputs be?
* what should the output be?
For a polygon clipping function, I assume that the inputs are two polygons,
A and B, and the output is a polygon C = A - (A intersect (Not B)). Am I
right? Does anyone know of an article or comme-il-faut way to do this?
Presumably, polygon clipping would depend on a more general polygon
intersection function, if this does not already exist?
Bjarke Christensen
------------------------------
Message: 21
Date: Mon, 22 Feb 2010 08:28:54 +0100 (CET)
From: Roger Bivand <Roger.Bivand at nhh.no>
To: Markus Loecher <markus.loecher at gmail.com>
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] clipping polygons, overlaying big polygon on
polygons
Message-ID: <alpine.LRH.2.00.1002220807280.32328 at reclus.nhh.no>
Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII
On Sun, 21 Feb 2010, Markus Loecher wrote:
I should have been much more specific with my question and supplied an
example.
In the meantime I believe I found the solution in the PBSmapping package,
this time I will provide an example that illustrates what I am trying to
achieve:
#Get the boundaries for Manhattan:
library(maps);
library(PBSmapping);
Manhattan <- map('county', 'new york,new york', plot=F);
#define a few points in Manhattan and compute the corresponding Voronoi
regions:
XY <- cbind.data.frame( PID = 1, POS = 1:10,
X =
c(-73.97220,-73.95456,-73.97906,-73.98685,-73.99968,-73.96716,-73.97127,-73.96524,-73.98322,-74.00587),
Y =
c(40.74983,40.77402,40.75022,40.73952,40.72842,40.75677,40.75112,40.75513,40.76271,40.74026)
); polys <- calcVoronoi(as.PolySet(XY, projection="LL") ); #plot the Voronoi regions: plotMap(polys) #It would be a much more meaningful plot to clip these Voronoi regions to the Manhattan boundaries: m <- as.PolySet (cbind.data.frame(PID = 1, POS = 1:length(Manhattan$x),
X =
Manhattan$x, Y = Manhattan$y)) plotMap(joinPolys(m, polys, "INT" ), xlim = Manhattan$range[1:2], ylim = Manhattan$range[3:4]); So the key operation I was looking for was simply the intersection
provided
by the function (joinPolys(..., "INT" ) from PBSmapping.
Yes, but sadly this does not solve the problem if you are going to use the output beyond hobby or education, as PBSmapping also uses the GPC code. Attempts were made to persuade the author to license his code in a more helpful way, but he has left his employing university, and they have retained the rights to his software. Although the PBSmapping maintainers claim to have permission to ship GPC under GPL, I would try to stay as far away from it as possible. Once code using it gets into the workflow, it is really hard to remove. The spatstat maintainers are working hard to make users aware of the problems, which will not go away, I'm afraid. If you'd like to try rgeos, please checkout or download from R-Forge. The gpc.poly-shadowing classes and methods now provide an intersection operation. Kicking the GPC habit is hard, but must be done, otherwise we risk spreading downstream dependencies. Before long CRAN may have colour-coded license dependencies, and organisations can already block the installation of CRAN packages that have non-free (as in speech) licenses or that depend upstream on non-free packages. The akima package is another worry, but has a much clearer license that does permit research under noncommercial use. GPC only permits "non-commercial use of GPC (for example: private / hobbyist / education)" while "commercial research" requires a commercial use license. Any contract research could fall under "commercial research", because it generates income. I would be very grateful for help with rgeos ... Roger