1

This question is similar to Create polygons from point dataset, where each polygon contains 3 points from the dataset, and in particular my problem is actually described in this comment to that question. But the problem were never solved there.

I have three sets of points. I want to partition the area which is covered by the points into polygons. Each polygon must have at least one point from each set of points.

The more polygons the better solution. I use R, so solutions written in R are preferred.

Ideally, the solution would be general enough to apply to a situation where there were two sets of points, or four sets of points.

My data is available here, and with this code, a map with the points can rendered:

library(ggmap)
load(url("http://hansekbrand.se/gis.RData")) ## this loads "gis" and "my.map"
my.map + geom_point(aes(colour = Generation, x = lon, y = lat),
                   size = 0.01, data = gis)

The three sets of points

The type of point is given by the variable Generation in the object gis

str(gis)
'data.frame':   875 obs. of  3 variables:
$ lon       : num  -1.52 -2.42 -2.35 -1.57 -1.53 ...
$ lat       : num  12.4 13.6 13.2 13.6 13.3 ...
$ Generation: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
Community
  • 1
Hans Ekbrand
  • 143
  • 5
  • You need to supply more criteria if you would like a unique solution. For instance, it's easy to produce solutions where the number of polygons is as large as logically possible (equal to the smallest number of points), but those polygons might need to be very tortuous. Do this by selecting one point from each set, connecting them in a tree that does not intersect any other point, and buffer that tree by such a small radius that the buffer contains no additional points. Remove the buffer from the region. Apply this procedure recursively to the remaining region and remaining points. – whuber Dec 10 '16 at 16:41
  • I do not need a unique solution, the ultimate use of the partitioning is to aggregate data for each polygon and use the polygons in a multi-level regression model. Each point represents a measurement, and what I want is to have comparable measurements in the same polygon, ideally the "most comparable ones", but as long as the points are close (in a loose sense), the partitioning will be useful. – Hans Ekbrand Dec 11 '16 at 12:37
  • You have just specified two more criteria that are absent in your question: "most comparable" and "close." Without those, literally any partitioning of the points into groups of three will work. To make progress, then, it's important that you quantify what you mean by "most comparable" and that you stipulate how you plan to trade off "comparability" for "closeness" when evaluating potential solutions. – whuber Dec 11 '16 at 16:05

1 Answers1

1

Just thinking out loud, not sure how to do it...

  1. Draw Voronoi diagram (package{deldir}), something like a picture below.
  2. Build polygons from obtained segments. (tile.list(){deldir})
  3. Create a field which records the variable Generation for each polygon.
  4. Define neighbors by poly2nb()
  5. Randomly select some neighbor pairs. (...how many?)
  6. Determine if neighbor pairs have same Generation variable.
  7. Merge neighbor pairs if they contain different Generation.
  8. Check Generation field of all polygons. Repeat steps 4 to 7 until all polygons have at least one point from every Generation.

voronoi

Hans Ekbrand
  • 143
  • 5
Kazuhito
  • 30,746
  • 5
  • 69
  • 149
  • 1
    I like your suggestions, but i found point 7 and 8 rather unclear, so I edited the wording. like this: "7 .Merge neighbor pairs if they contain different Generation. 8. Check Generation field of all polygons. Repeat steps 4 to 7 until all polygons have at least one point from every Generation." Please tell if that's what you meant. (The edit might will not visible until it is peer reviewed). – Hans Ekbrand Dec 10 '16 at 13:25
  • @HansEkbrand Yes, it's now clear and much better. Glad my rather poor idea went into good hands. – Kazuhito Dec 10 '16 at 13:53