Skip to content

R-sig-Geo Digest, Vol 43, Issue 14

13 messages · Nicholas Lewin-Koh, White.Denis@epamail.epa.gov, Roger Bivand +4 more

Original
Nicholas Lewin-Koh
# 
Hi Tim,
I am not quite sure what you are getting at here. Do you want to
intersect
polygons and then select the set of lines that form the outer perimeter?
Do you wan the convex hull of a set of polygons. I guess I have been out
of the
GIS world to long. It seems to me that this would be something easy to
solve,
just tedious iteration of the polygon coordinates and some
triangulation.

Nicholas
Tim Keitt
# 
Hi Nic,

The convex hull would be fast and easy to compute (there's existing
code in R). I want the ordinary hull which is the set of arcs forming
the perimeters (inside and out). My crude and very slow solution was
to convert all the polygons (in this case hexagons on a lattice) into
their constituent arcs and then for each arc count how many times it
occurs in the set (requires slightly fuzzy matching of points). Arcs
that occur more than once are removed. The remaining arcs form the
hull. Runs in about 20 minutes with a  few hundred hexagons.
Sufficient for the moment.

THK
On 3/16/07, Nicholas Lewin-Koh <nikko at hailmail.net> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Nicholas Lewin-Koh
    

    

      
      #
    

  


  

      
Hi Tim,
You could compute the convex hull first, and then iterate from points
on the convex hull. That should be much faster already, especially since
hexagons are convex and the perimeter will be locally convex around all
the
points touching the convex hull. You could do a variation
of the "monotone pieces" algorithm that is used in computational
geometry.
But this is a simpler problem. Are there cases with interior holes?

I have been meaning to write something like this for hexbin for a while.
There
are many cases where it would be nice to find approximations to the
density contours
and a quick and dirty way is to threshold the hexagon counts, find the
hull and
smooth the perimeter.

Nicholas

On Fri, 16 Mar 2007 08:34:20 -0500, "Tim Keitt" <tkeitt at gmail.com> said:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Tim Keitt
    

    

      
      #
    

  


  

      
Ya, there can be pretty arbitrary holes. The situation is a large
hexagonal lattice where some fraction of cells have q-values < 0.05.
Now draw lines around those cells following the hexagonal edges, but
omitting any interior arcs. I can send you an example solution. The
images are bit big for the list.

THK

        
On 3/16/07, Nicholas Lewin-Koh <nikko at hailmail.net> wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      White.Denis@epamail.epa.gov
    

    

      
      #
    

  


  

      
The thread below and many others in R-sig-Geo raise questions about
future directions.  In reinventing GIS there are a whole list of
capabilities and functions that would be helpful.  Some that I have
noticed include,

Topological representation to enable
      Planar enforcement of boundary integrity of polygon tessellations
      "Dissolving" interior edges easily as in the thread below
Large problem computational geometry functions
      Identify many points inside of many polygons
      Intersections/overlays of two sets of many polygons
      Distances between all pairs of many polygons

Are there members of the R Geo community working on any of these?
Are these issues seen as an exclusive focus of commercial GIS?
Are there discussions about these issues at relevant conferences?

(I will be at AAG in San Francisco and would be happy to meet with
others if there is interest.)

r-sig-geo-bounces at stat.math.ethz.ch wrote on 2007-03-16 07:52:59:

            

      
      
      
    

        
since

            

      
      
      
    

        
all

            

      
      
      
    

        
while.

            

      
      
      
    

        
said:

            

      
      
      
    

        
forming

            

      
      
      
    

        
into

            

      
      
      
    

        
perimeter?

            

      
      
      
    

        
been out

            

      
      
      
    

        
easy to

            

      
      
      
    

        
<6262c54c0703150849qe60ab14nfef1eb3bf73dfb5d at mail.gmail.com>

            

      
      
      
    

        
and

            

      
      
      
    

        
the

            

      
      
      
    

        
arcs

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Nicholas Lewin-Koh
    

    

      
      #
    

  


  

      
Hi Tim,
Now I think I get it. What about the following:
1) threshold based on your qvalues
2) find all disconnected sets 
3) For each disconnected set 
  a) find the convex hull of the hexagon centers, this should be
  equivalent to the
     the convex hull of the hexagons themselves (convexity of the
     hexagons and a regular lattice)
  b) iterate from each (center)point on the hull, using only points in
  the set and test for an exterior neighbor,
     you have to be careful here due to holes, just keep a binary matrix
     of each hexagons exterior faces.
  c) invert the polygon and find the boundaries of interior holes
4) put it all together saving just the vertices of the hexagons sides on
the exterior

So you have already reduced your constant by 1/6 using the centers. I
think this should be no worse
than the convex hull itself times some function of n representing the
fraction of exterior polygon
centers.

Now you already have something that works so this is just thinking about
if you have to do it again or
do some simulations.

Cheers
Nicholas
  

On Fri, 16 Mar 2007 10:55:06 -0500, "Tim Keitt" <tkeitt at gmail.com> said:

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 16 Mar 2007 White.Denis at epamail.epa.gov wrote:

        

            

      
      
      
    

        
One possibility was mentioned in a question by Hisaji Ono last year, and 
which could be done, is to build rgdal against a GDAL/OGR with GEOS 
support, and use the OGR functions such as:

virtual OGRGeometry * 	Union (const OGRGeometry *) const

and other similar ones. This would also involve trying to make a faster 
conversion between the OGR and internal sp representations, and may 
precipitate a redefinition of polygons in sp to be fully OGR compatible.

The sp S4 classes would need much better integration with OGR for this to 
work adequately on larger data sets, but we could profile it after 
prototyping.

            

      
      
      
    

        
Thanks for taking this initiative! As one of the guilty parties, and also
at the AAG, count me in, at least for the prototyping - I think GEOS is a
realistic place to start because it is being used by others and maintained
actively.

Having tried things out, we may find that this is too much GIS, and that 
either using a GIS or a database interface like PostGIS or aRT/Terralib is 
a better idea. But at least we'd know.

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Tim Keitt
    

    

      
      #
    

  


  

      
Denis,

I was just about to reply to Nicholas that my ideal solution would be
to build a graph data structure linking adjacent elements and find the
the hull that way -- I think you saw the solution before I did. (I'm
starting to grok that all of this is already well known in the vector
GIS world where topology has been thoroughly studied.)

I'd be happy to discuss with folks the idea of using the Boost Graph
Library to implement topological relationships for spatial objects.
We've already produced code that links the BGL to GDAL -- specify an
adjacency rule and then run eg cost path analysis on GDAL rasters
(look for this on sourceforge soon). It would also be possible to
implement something similar with vector objects.

My OSGeo summer of code idea is to link the BGL to postgis. I wasn't
originally thinking about topology per se, but it would achieve much
the same thing.

If people have ideas about how to fund these sorts of things, I'd be
happy to hear about it. The implementation end is rather far from my
core research. Its mostly been side projects.

THK

        
On 3/16/07, White.Denis at epamail.epa.gov <White.Denis at epamail.epa.gov> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      elijah wright
    

    

      
      #
    

  


  

          

      
      
      
    

        
I've thought a bit about this project, too, as it butts up against some of 
my own work.

My experience with graph libraries has not been altogether positive.  ;) 
They are rather difficult to coerce into handling very large data sets.

I think it would be handy to have a layer to pump graph data through BGL 
into a postgis store.  It would be even better if that postgis store 
looked something similar to the store used as a backend by the 
Redland/Raptor xml graph serialization bits.  :-)  Just a suggestion for a 
place to start.... I know the redland/raptor bits are extremely 
scale-friendly for attributed graph data.

--elijah

  
  


  


      

    

    

      
      

        


  

        

      Tim Keitt
    

    

      
      #
    

  


  

      
On 3/16/07, elw at stderr.org <elw at stderr.org> wrote:

            

      
      
      
    

        
Our approach is not to store the vertex/edge information in the graph
library, but rather to model it with iterator interfaces. Generic
programming methods make this possible. Essentially, GDAL becomes the
data storage mechanism and the graph is implicit. The limits in this
case are likely the algorithmic complexity rather than data storage
(the step of building a graph data structure is entirely avoided). Of
course there may be better alternatives that are custom built for the
problem.

Vertex/edge iterators could also be backed by postgis or any other
data source with sufficient structure to specify adjacency. The
enormous value from my perspective is that effort coding algorithms
becomes reusable in many contexts.

THK

            

      
      
      
    

        
I've never heard of it, but like I say, this is pretty far from my core areas.

THK

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Edzer Pebesma
    

    

      
      #
    

  


  

      
Although I'm interested many of the issues below, and also feel the 
attraction of the develop-it-ourself, my personal agenda is a bit more 
modest. Once we start with GEOS support (an idea mentioned every 6 
months on this list), the next thing we need for dealing with huge or 
massive data sets is a clever indexing structure. After that, full 
topology. I find it hard to imagine where the resources should come from 
for all these things that are available in OSGeo next door -- we can't 
settle with proof of concepts, but want quality stuff. OTOH, one day a 
student may come in who just finished two of these projects.

My primary goal is to get sp stable, and adopted by more of the numerous 
packages for analyzing spatial data in R. And yes, I did introduce an 
instabillity this week -- try [ on a SpatialPolygonsDataFrame with NA 
values in the row index. There's a lot of code, much of which is little 
tested and/or not very clean. Another goals is to get a good and smooth 
system where R works as a back-end in interopable systems, possibly 
using PostGIS to transfer data. Finally: further proof that R is a 
increasingly wonderful system for analyzinig spatial data.
--
Edzer

        
White.Denis at epamail.epa.gov wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Gabor Grothendieck
    

    

      
      #
    

  


  

      
In my opinion the first item on the agenda should be to
get good documentation in place including many examples.

        
On 3/16/07, Edzer J. Pebesma <e.pebesma at geo.uu.nl> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    
1 day later

    

      
      

        


  

        

      Tim Keitt
    

    

      
      #
    

  


  

      
I agree GEOS support would be welcome. I'd be happy to help out if
anyone takes this on.

I've been thinking of a possibly quicker (devel time) solution lately
since I've been using PostGIS a lot. One could have a postgis package
that copies 'sp' object into postgis and then provides postgis
function equivalents that operate on the R proxy objects.

pgobj1 <- new('PostGISGeom', my.spatial.polygons)
pgobj2 <- somePostGISOp(pgobj1, ...)

res.spatial.polygons <- PGObj2SPObj(pgobj2)  # or 'as' generic

The pgobj? objects could be first class 'sp' members if it is possible
to subclass the 'sp' geometry classes -- may not be possible as I
don't think you can overload the S4 slot call.

THK

        
On 3/16/07, Edzer J. Pebesma <e.pebesma at geo.uu.nl> wrote: