Message Boards Message Boards

GROUPS:

Voronoi and borders

Posted 6 years ago
9199 Views
|
5 Replies
|
5 Total Likes
|

Hello I create a voronoi region from points that locate the cities of a country. I also creates the border of this country. I want to erase the parties to External voronoi polygon borders. What is the solution? I try RegionIntersection and RegionDifference, but this don't work. Thank you in advance. My little program :

André Dauphiné

   ClearAll["Global`*"]

   pays = "France";

   ny = DialogInput[
      DynamicModule[{name = ""}, Column[{"How much cities?",
         InputField[Dynamic[name], String],
   ChoiceButtons[{DialogReturn[name], DialogReturn[]}]}]]];

   ny = ToExpression[ny];

   coord = Take[Table[Reverse[CityData[c, "Coordinates"]], {c, CityData[{All, pays}]}], ny];

   border = RegionBoundary[ConvexHullMesh[coord]]

   voron = VoronoiMesh[coord]
5 Replies

Unfortunately I do not have much experience in working with regions. But I find your question so interesting that I am trying now a "solution":

ClearAll["Global`*"]

pays = "France";
ny = ToExpression@
   DialogInput[
    DynamicModule[{name = ""}, 
     Column[{"How much cities?", InputField[Dynamic[name], String], 
       ChoiceButtons[{DialogReturn[name], DialogReturn[]}]}]]];
coord = Take[
   Table[Reverse[CityData[c, "Coordinates"]], {c, 
     CityData[{All, pays}]}], ny];
border = ConvexHullMesh[coord];   (* !!! *)
voron = VoronoiMesh[coord];
chm = ConvexHullMesh @@@ MeshPrimitives[voron, 2];

(* here is my idea: the individual intersection of the complete
region with every Voronoi tile: *)
ri = RegionIntersection[border, #] & /@ chm;
grLines = MeshPrimitives[#, 1] & /@ ri;
gr = Graphics[grLines]
DiscretizeGraphics[gr]

While this is quite a crude approach, I am sure there are much more elegant solutions. Maybe someone is going to give a hint.

Cheers Henrik

Thanks, It's OK André

Hi Andre,

it is not OK!

The solution in my response was just given in terms of your question. But a general remark must be made (as I found during a sleepless night): A Voronoi diagram basically deals with distances. Angular coordinates are therefore not appropriate here. If one wants to do it right, one has to perform a transformation into grid coordinates first. I tried, and my little code now looks like this:

ClearAll["Global`*"]

(* just 50 french cities are taken into account: *)
coord = Take[
   Table[CityData[c, "Coordinates"], {c, CityData[{All, "France"}]}], 
   50];

(* transformation from angular coordinates to grid-coordinates: *)
coordxy = 
  GeoGridPosition[
      GeoPosition[#, "WGS84"], {"UTMZone31", 
       "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}][[
     1]] & /@ coord;

(* now the voronoi tessellation is made: *)
border = ConvexHullMesh[coordxy];
voron = VoronoiMesh[coordxy];
chm = ConvexHullMesh @@@ MeshPrimitives[voron, 2];
ri = RegionIntersection[border, #] & /@ chm;
grLines = MeshPrimitives[#, 1] & /@ ri;

(* procedure for transformation graphics lines to GeoPath: *)
makeGeoPath[lineList_] := Module[{pts, geoPts},
  pts = lineList[[;; , 1, 1]];
  geoPts = 
   GeoPosition[GeoGridPosition[#, "UTMZone31"]] & /@ 
    Append[pts, pts[[1]]];
  GeoPath[geoPts]
  ]

(* show the result: *)
GeoGraphics[{Blue, makeGeoPath /@ grLines, Red,
   GeoMarker[GeoPosition /@ coord, \[FilledCircle]]}, 
 ImageSize -> Full]

The output then is: enter image description here

With using simply the angular coordinates, the result is just wrong: enter image description here

Cheers Henrik

Well, I would like to add a final remark: Instead of the convex hull of the cities point cloud the real border polygon can/should be used. Then the whole thing migth even make some sense. My code:

ClearAll["Global`*"]

(* coordinates of 200 French cities: *)
cityCoords = (CityData[#, "Coordinates"] & /@ 
     CityData[{All, "France"}])[[;; 200]];
frenchPolygon = CountryData["France", "Polygon"];
borderCoords = frenchPolygon[[1, 1]];

(* transformation from angular coordinates to grid-coordinates: *)
geoGridPos = 
  GeoGridPosition[
     GeoPosition[#, "WGS84"], {"UTMZone31", 
      "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}][[
    1]] &;

cityCoordsxy = geoGridPos /@ cityCoords;
borderCoordsxy = Map[geoGridPos, borderCoords, {2}];

(* now the voronoi tessellation is made: *)
border = BoundaryDiscretizeGraphics[Graphics[Polygon[borderCoordsxy]]];
voron = VoronoiMesh[cityCoordsxy];
chm = ConvexHullMesh @@@ MeshPrimitives[voron, 2];
ri = RegionIntersection[border, #] & /@ chm;
grLines = MeshPrimitives[#, 1] & /@ ri;

(* procedure for transformation graphics lines to GeoPath: *)
makeGeoPath[lineList_] := Module[{pts, geoPts},
  pts = lineList[[;; , 1, 1]];
  geoPts = 
   GeoPosition[GeoGridPosition[#, "UTMZone31"]] & /@ 
    Append[pts, pts[[1]]];
  GeoPath[geoPts]
  ]

(* show the result: *)
GeoGraphics[{Green, frenchPolygon, Red, 
  GeoMarker[GeoPosition /@ cityCoords, \[FilledCircle]], Blue, 
  makeGeoPath /@ grLines}, ImageSize -> Full]

enter image description here

MerryChristmas /@ wolframCommunity

Splendid Can I have your adress mail

Bonnes et joyeuses fetes

André Dauphiné

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract