# Voronoi and borders

Posted 7 years ago
9390 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
Sort By:
Posted 7 years ago
 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 7 years ago
 Thanks, It's OK André
Posted 7 years ago
 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: With using simply the angular coordinates, the result is just wrong: 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] MerryChristmas /@ wolframCommunity