Message Boards Message Boards

2
|
11692 Views
|
5 Replies
|
5 Total Likes
View groups...
Share
Share this post:

Voronoi and borders

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]
POSTED BY: André Dauphiné
5 Replies

Splendid Can I have your adress mail

Bonnes et joyeuses fetes

André Dauphiné

POSTED BY: André Dauphiné

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

POSTED BY: Henrik Schachner

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

POSTED BY: Henrik Schachner

Thanks, It's OK André

POSTED BY: André Dauphiné

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

POSTED BY: Henrik Schachner
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