# Plotting Voronoi diagram for a country

Posted 1 year ago
2952 Views
|
25 Replies
|
30 Total Likes
|
 Dear all,I have geographical coordinates and rainfall values for 151 points (stations) over Turkey. How do I plot the Voronoi diagram over this country with Mathematica? Sations = {"ADANABOLGE", "ADIYAMAN", "AFYONKARAHİSARBÖLGE", "AĞRI", "AHLAT", "AKÇAABAT", "AKÇAKOCA", "AKSARAY", "AKŞEHİR", "AMASRA", "AMASYA", "ANKARABÖLGE", "ARDAHAN", "ARTVİN", "AYVALIK", "BAFRA", "BALIKESİRGÖNEN", "BANDIRMA", "BARTIN", "BAŞKALE", "BAYBURT", "BERGAMA", "BEYPAZARI", "BEYŞEHİR", "BİLECİK", "BİNGÖL", "BOĞAZLIYAN", "BOLU", "BOLVADİN", "BOZÜYÜK", "BURDUR", "BURSA", "ÇANAKKALE", "ÇANKIRI", "ÇEMİŞGEZEK", "CEYHAN", "CİHANBEYLİ", "ÇORLU", "ÇORUM", "DENİZLİ", "DEVELİ", "DİNAR", "DİVRİĞİ", "DİYARBAKIRHAVALİMANI", "DÖRTYOL", "DURSUNBEY", "DÜZCE", "EDİRNE", "EĞİRDİR", "ELAZIĞBÖLGE", "ELMALI", "EMİRDAĞ", "ERCİŞ", "ERDEMLİ", "EREĞLİ", "ERZİNCAN", "ERZURUMHAVALİMANI", "FLORYA", "GAZİANTEP", "GEMEREK", "GEYVE", "GİRESUN", "GÖKÇEADA", "GÖKSUN", "GÜMÜŞHANE", "GÜNEY", "HADİM", "HINIS", "HOPA", "HORASAN", "IĞDIR", "ILGIN", "İNEBOLU", "İSKENDERUN", "ISPARTA", "İSPİR", "İZMİRBÖLGE", "KAHRAMANMARAŞ", "KARAİSALI", "KARAMAN", "KARAPINAR", "KARS", "KASTAMONU", "KASTAMONUBOZKURT", "KAYSERİBÖLGE", "KEBAN", "KELES", "KİLİS", "KIRKLARELİ", "KIRŞEHİR", "KIZILCAHAMAM", "KOCAELİ", "KONYAHAVALİMANI", "KORKUTELİ", "KOZAN", "KULU", "KÜTAHYA", "MALATYA", "MALAZGİRT", "MANAVGAT", "MERZİFON", "MUĞLA", "MUŞ", "NALLIHAN", "NEVŞEHİR", "NİĞDE", "ORDU", "ÖZALP", "PALU", "POLATLI", "RİZE", "RİZEPAZAR", "SAKARYA", "SALİHLİ", "SAMANDAĞ", "SAMSUNBÖLGE", "SARIKAMIŞ", "SARIYER", "SARIYERKUMKÖYKİLYOS", "SARIZ", "ŞEBİNKARAHİSAR", "SENİRKENT", "SEYDİŞEHİR", "SİİRT", "ŞİLE", "SİMAV", "SİNOP", "SİVAS", "SIVRIHISAR", "SOLHAN", "TAVŞANLI", "TEFENNİ", "TEKİRDAĞ", "TERCAN", "TOKAT", "TORTUM", "TOSYA", "TUNCELİ", "ULUBORLU", "ULUDAĞ", "UŞAK", "UZUNKÖPRÜ", "VANBÖLGE", "YALOVA", "YATAĞAN", "YOZGAT", "YUMURTALIK", "YUNAK", "ZARA", "ZİLE", "ZONGULDAK"}; lat = {37., 37.75, 38.75, 39.72, 38.77, 41.02, 41.08, 38.38, 38.35, 41.75, 40.65, 39.95, 41.12, 41.18, 39.32, 41.57, 40.1, 40.35, 41.63, 38.05, 40.25, 39.12, 40.17, 37.68, 40.15, 38.88, 39.2, 40.73, 38.72, 39.92, 37.72, 40.18, 40.15, 40.6, 39.07, 37.03, 38.65, 41.17, 40.55, 37.78, 38.38, 38.07, 39.37, 37.9, 36.85, 39.58, 40.83, 41.67, 37.87, 38.67, 36.75, 39.02, 39.03, 36.62, 37.5, 39.75, 39.92, 40.98, 37.07, 39.18, 40.52, 40.92, 40.2, 38.02, 40.47, 38.15, 36.98, 39.37, 41.4, 40.05, 39.92, 38.28, 41.98, 36.58, 37.77, 40.48, 38.43, 37.6, 37.27, 37.18, 37.72, 40.62, 41.37, 41.95, 38.73, 38.8, 39.92, 36.72, 41.73, 39.15, 40.47, 40.78, 37.87, 36.75, 37.45, 39.1, 39.42, 38.35, 39.15, 36.78, 40.87, 37.22, 38.73, 40.18, 38.58, 37.97, 40.98, 38.67, 38.72, 39.58, 41.03, 41.18, 40.78, 38.48, 36.08, 41.28, 40.33, 41.17, 41.25, 38.48, 40.3, 38.1, 37.42, 37.92, 41.18, 39.08, 42.02, 39.75, 39.45, 38.97, 39.55, 37.32, 40.98, 39.78, 40.3, 40.3, 41.02, 39.12, 38.08, 40.13, 38.68, 41.27, 38.5, 40.65, 37.35, 39.82, 36.77, 38.82, 39.9, 40.3, 41.45}; lon = {35.33, 38.28, 30.53, 43.05, 42.5, 39.58, 31.13, 34.08, 31.42, 32.38, 35.83, 32.88, 42.72, 41.82, 26.7, 35.92, 27.65, 27.97, 32.33, 44.02, 40.23, 27.18, 31.92, 31.72, 29.98, 40.48, 35.25, 31.52, 31.05, 30.03, 30.28, 29.07, 26.42, 33.62, 38.92, 35.82, 32.93, 27.8, 34.95, 29.08, 35.5, 30.17, 38.12, 40.23, 36.22, 28.63, 31.17, 26.57, 30.83, 39.23, 29.92, 31.15, 43.35, 34.3, 34.05, 39.5, 41.27, 28.75, 37.38, 36.07, 30.3, 38.4, 25.9, 36.5, 39.47, 29.07, 32.47, 41.7, 41.43, 42.17, 44.05, 31.92, 33.77, 36.17, 30.55, 41., 27.17, 36.93, 35.07, 33.22, 33.55, 43.1, 33.78, 34.02, 35.48, 38.75, 29.07, 37.12, 27.23, 34.17, 32.65, 29.93, 32.48, 30.2, 35.82, 33., 29.97, 38.32, 42.53, 31.43, 35.33, 28.37, 41.48, 31.35, 34.67, 34.68, 37.9, 43.98, 39.97, 32.15, 40.52, 40.88, 30.42, 28.13, 35.97, 36.3, 42.57, 29.05, 29.03, 36.5, 38.42, 30.55, 31.83, 41.95, 29.61, 28.98, 35.17, 37.02, 31.53, 41.07, 29.5, 29.77, 27.55, 40.38, 36.57, 41.55, 34.03, 39.55, 30.45, 29.08, 29.4, 26.68, 43.38, 29.27, 28.13, 34.8, 35.78, 31.73, 37.75, 35.75, 31.8}; rain = {65.10364583333333, 69.24427083333333, 35.316875, 44.31385416666667, 48.824375, 60.07927083333334, 90.26499999999999, 30.4459375, 47.152499999999996, 84.07041666666667, 39.665104166666666, 33.730937499999996, 46.18989583333334, 59.59437500000001, 65.39645833333333, 66.26520833333333, 56.8996875, 60.803437499999994, 86.75604166666668, 37.60124999999999, 37.283229166666665, 62.348958333333336, 34.07760416666667, 42.45520833333334, 37.88489583333333, 83.88291666666665, 33.40291666666666, 46.26270833333333, 33.0103125, 40.45520833333333, 35.920833333333334, 59.19114583333333, 54.93541666666667, 33.924166666666665, 54.429999999999986, 68.190625, 28.725729166666664, 48.099687499999995, 38.19916666666666, 50.333125, 33.5215625, 38.09729166666667, 33.66562499999999, 47.78999999999999, 80.06635416666668, 47.6740625, 68.9203125, 49.01083333333332, 70.82197916666667, 36.59572916666666, 42.33979166666666, 34.72958333333332, 37.396875, 61.012187499999996, 27.422604166666662, 31.66, 34.3028125, 54.60708333333333, 59.14447916666666, 36.41260416666667, 53.246249999999996, 104.93041666666666, 64.91343749999999, 54.14729166666666, 38.96833333333334, 48.09666666666667, 58.42135416666667, 50.02572916666667, 186.38958333333332, 34.09375, 21.917500000000004, 36.49583333333333, 84.91864583333333, 63.12968749999999, 43.601875, 39.66906250000001, 74.27645833333334, 75.86124999999998, 78.48114583333333, 30.548333333333332, 27.080729166666668, 40.853541666666665, 41.176770833333336, 102.25729166666665, 34.189791666666665, 33.369791666666664, 62.652604166666656, 53.38614583333333, 46.46583333333333, 33.14489583333333, 47.733333333333334, 68.51447916666666, 27.64041666666667, 33.63822916666667, 72.35229166666667, 33.57364583333333, 45.55354166666667, 34.13260416666666, 39.54677083333333, 129.77604166666666, 37.36052083333333, 109.60208333333331, 67.33604166666666, 28.4934375, 36.56010416666666, 28.66677083333333, 86.56166666666665, 41.82364583333333, 50.02427083333333, 30.500416666666666, 186.503125, 170.72708333333335, 70.54874999999998, 46.05270833333334, 89.77739583333333, 57.81166666666666, 51.84885416666666, 70.90697916666667, 70.74885416666667, 46.635833333333345, 50.375104166666674, 56.90666666666666, 63.61427083333333, 65.45572916666667, 73.02739583333334, 69.20229166666667, 56.49010416666667, 38.18552083333333, 33.69833333333334, 59.517916666666665, 40.5925, 39.89760416666666, 49.4359375, 37.839166666666664, 37.82166666666667, 40.09354166666666, 39.73760416666667, 75.17218749999999, 52.24187499999999, 123.11343749999999, 47.08375, 55.36291666666667, 34.23145833333333, 62.806562500000005, 62.40125, 52.1103125, 80.3328125, 38.78125, 45.807291666666664, 39.19291666666667, 102.12302083333334}; latLngs = Transpose[{lat, lon}]; toMercator[{lat_, lng_}] := {lng, Log[Abs[Sec[lat*Degree] + Tan[lat*Degree]]]/Degree}; mercPoints = toMercator /@ latLngs; Show[CountryData["Turkey", {"Shape", "Mercator"}], Frame -> True, Epilog -> {PointSize[0.01], Red, Point[mercPoints]}] 
Answer
25 Replies
Sort By:
Posted 1 year ago
 Hi Mohammad,There are two ; missing in the definition of makeGeoPath. makeGeoPath[lineList_] := Module[{pts, geoPts}, pts = lineList[[;; , 1, 1]]; geoPts = GeoPosition[GeoGridPosition[#, "UTMZone31"]] & /@ Append[pts, pts[[1]]]; GeoPath[geoPts]] 
Answer
Posted 1 year ago
 Hello Muhammad,here comes one more attempt. Let stations be the list of the location names, etc., then: turkeyPolygon = CountryData["Turkey", "Polygon"]; borderCoords = turkeyPolygon[[1, 1]]; geoGridPos = GeoGridPosition[ GeoPosition[#, "WGS84"], {"UTMZone36", "CentralScaleFactor" -> 0.9996, "GridOrigin" -> {500000, 0}}][[1]] &; cityCoordsxy = geoGridPos /@ cityCoords; borderCoordsxy = Map[geoGridPos, borderCoords, {2}]; border = BoundaryDiscretizeGraphics[Graphics[Polygon[borderCoordsxy]]]; voron = VoronoiMesh[cityCoordsxy]; chm = BoundaryMeshRegion /@ MeshPrimitives[voron, 2]; colors = ColorData["AvocadoColors"] /@ Rescale[rain]; cell0 = Thread[{stations, cityCoordsxy, rain, colors}]; cell1 = {#1, #2, #3, #4, First@Select[chm, Function[{poly}, #2 \[Element] poly]]} & @@@ cell0; cells = <|"station" -> #1, "coordxy" -> #2, "rain" -> #3, "color" -> #4, "polymesh" -> RegionIntersection[border, #5]|> & @@@ cell1; makeGeoCell[cell_] := Module[{pts}, pts = First /@ MeshPrimitives[cell["polymesh"], 0]; {cell["color"], Polygon[pts], PointSize[Medium], Red, Point@cell["coordxy"]}] GraphicsRow[{Graphics[makeGeoCell /@ cells], BarLegend[{"AvocadoColors", MinMax[rain]}]}, ImageSize -> 800] 
Answer
Posted 1 year ago
 In the result of my approach above I am using just Graphics and not GeoGraphics, because I have no idea on how to draw a Polygon within the GeoGraphics-environment. Therefore I regard this as sub-optimal! We have GeoRegionValuePlot (which would be ideal for this task) but there is no GeoRegion I could convert a Polygon into - I am talking about an opaque Polygon with a definite color. Or at least I could not find out how to do it. What am I missing?Apart from this one can play around with the given data, and there are - in particular in comparison with the lengthy code above - very easy and nice ways to present them, e.g.: wd = WeightedData[GeoPosition[cityCoords], rain]; GeoHistogram[wd, ImageSize -> 700, PlotLegends -> Automatic] GeoSmoothHistogram[wd, ImageSize -> 700] GeoBubbleChart[wd, ImageSize -> 700] @Alex Teymouri : Thank you for your kind words!
Answer
Posted 11 months ago
 ... too late I found out that ListDensityPlot[..., InterpolationOrder -> 0,...] is basically doing the same job (in a "quick and dirty way") - as a one-liner: ListDensityPlot[MapThread[Append, {cityCoordsxy, rain}], InterpolationOrder -> 0, ColorFunction -> "AvocadoColors", AspectRatio -> Automatic, ImageSize -> 700, PlotRange -> Flatten[{CoordinateBounds@border, All}, 1], Frame -> None, Epilog -> {Red, Point[cityCoordsxy]}, RegionFunction -> ({#1, #2} \[Element] border &), PlotLegends -> Automatic] 
Answer
Posted 1 year ago
 You have a similar question by André Dauphiné and a response by Henrik Schancher A. Dauphiné
Answer
Posted 1 year ago
 @André Dauphiné: Thanks for the reference and greetings!The respective discussion can be found here:https://community.wolfram.com/groups/-/m/t/410510
Answer
Posted 9 months ago
 Dear Muhammad,I am afraid you have to be much more specific on what is wanted - I can only guess that you want a connection from any station with "Strength" s to all other stations with "Strength" #1, "Name" -> #2] &, {Strength, Stations}]; comps0 = Tuples[stats[[ ;; ]], 2]; comps1 = Select[comps0, First[#]["Strength"] > Last[#]["Strength"] &]; grRules = DirectedEdge[First[#]["Name"], Last[#]["Name"]] & /@ comps1; Graph[grRules, VertexLabels -> Automatic] The resulting Graph then has 11325 edges - plotting this is quite meaningless. Regards -- Henrik
Answer
Posted 11 months ago
 Hello Mohammad,the area of the polygons can be calculated quite easily. I am using the content of cells from my code in this post, and make then a Dataset from it and extend it with the respective "polyArea". cellDS = Dataset[cells][All, <|"station" -> #station, "coordxy" -> #coordxy, "rain" -> #rain, "color" -> #color, "polymesh" -> #polymesh, "polyArea" -> Quantity[Area[#polymesh], "Meters"^2]|> &] Because the coordinates of the polygons are already in [meters] (thanks to GeoGridPosition), the result comes directly in [meters^2]. Of course this is just an approximation, since the calculation is done for plane geometry, but it should be tolerable. You can convince yourself simply like so: (* true area: *) Entity["Country", "Turkey"]["Area"] (* approximated area (sum of all polygons): *) Total@cellDS[All, "polyArea"] Regards -- Henrik
Answer
Posted 11 months ago
 Ok, then you did not evaluate eventually a line like stations = ... Anyway - please find attached my full notebook. Regards -- Henrik Attachments:
Answer
Posted 9 months ago
 Alex,Looks like those examples use a subset of locations and position the vertices according to geo position. Using Henrik's answer make this change (* lat/lon reversed to get x/y position *) g = Graph[grRules, VertexLabels -> Automatic, VertexCoordinates -> Transpose[{lon, lat}]] Select some subset of stations Subgraph[g, {"ZONGULDAK", "YUNAK", "İSKENDERUN", "BALIKESİRGÖNEN", "ÇORUM", "MUĞLA", "KONYAHAVALİMANI", "IĞDIR"}, VertexSize -> Tiny, EdgeStyle -> Arrowheads[.015]] 
Answer
Posted 9 months ago
 Alex,How is "linear correlation coefficient" defined based on strength values?
Answer
Posted 9 months ago
 Hi Alex,Change the filter criteria compsStrong = Select[comps0, First[#]["Strength"] > 11.37 && Last[#]["Strength"] > 11.37 && First[#]["Name"] != Last[#]["Name"] &]; grRulesStrong = DirectedEdge[First[#]["Name"], Last[#]["Name"]] & /@ compsStrong Since all stations are not included, need an association of station to coordinate and use it to get the coordinates for the vertices in the graph stationCoordinates = AssociationThread[Stations, Transpose[{lon, lat}]]; vertexCoordinates = stationCoordinates /@ grRulesStrong[[All, 1]] // DeleteDuplicates Add color to vertices based on strength stationStrengths = AssociationThread[Stations, Strength]; vertexStrengths = (# -> stationStrengths[#] &) /@ (grRulesStrong[[All, 1]] // DeleteDuplicates); vertexColors = Thread[Rule[Keys@vertexStrengths, ColorData["TemperatureMap"] /@ Rescale[Values@vertexStrengths, MinMax@Values@vertexStrengths, {0.5, 1.0}]]] Construct graph gStrong = Graph[grRulesStrong, VertexLabels -> Automatic, VertexCoordinates -> vertexCoordinates, VertexStyle -> vertexColors]; It still has a large number of vertices and edges Through[{VertexCount, EdgeCount}[gStrong]] (* {72, 5112} *) Pick a random sample of stations to plot SeedRandom[314]; Subgraph[gStrong, VertexList@gStrong // RandomChoice[#, 10] &, VertexSize -> {"Scaled", .02}, EdgeStyle -> Arrowheads[.01], ImageSize -> 600] 
Answer
Posted 1 year ago
 Thank you so much, Henrik and André.I did step by step but for the last row of the code, I got an error. Please see the enclosed notebook. Additionally, may we have a different color for each region and a legend like the attached figure?I really appreciate your time. Attachments:
Answer
Posted 1 year ago
 Dear Henrik and Rohit,Awesome solution for a difficult question. It’s so helpful. Thank you :)
Answer
Posted 1 year ago
 Rohit and Henrik, You did a perfect and nice solution.I deeply appreciate your great help.
Answer
Posted 11 months ago
 Thanks, Henrik. That is a perfect suggestion.
Answer
Posted 11 months ago
 Hi Rohit and Henrik,I have an important question, but I don't know if it can be calculated or not.Is it possible to estimate the area of each polygon based on the total area of the country?Regards,
Answer
Posted 11 months ago
 Dear Henrik,That's a nice solution. I really appreciate your time :)Is it possible to add station names to the output table?
Answer
Posted 11 months ago
 I do not quite understand! Just call cellDS - and you see everything together in one (Dataset)-table! Why don't you just try it?! Extraction of station name and area into a simple list can be done with one simple command: Normal@Values@cellDS[All, {"station", "polyArea"}] 
Answer
Posted 11 months ago
 Dear Henrik, I got this kind of output:
Answer
Posted 9 months ago
 Hi Henrik and Rohit,If you remember, I was studying about rain gauge stations in Turkey. You and Rohit helped me a lot.Now for drawing a network for a few points, I found the below link: https://mathematica.stackexchange.com/questions/170612/create-a-network-using-dataFor 151 raingauge stations we have a data list like this: DataList={(lat,lon,elevation,m,tau,strength,rain)1, (lat,lon,elevation,m,tau,strength,rain)2, (lat,lon,elevation,m,tau,strength,rain)3,…, (lat,lon,elevation,m,tau,strength,rain)151} Is it possbile to extract a network for this data? The data is attached.I never forget your kindness.RegardsMohammad Attachments:
Answer
Posted 9 months ago
 Thank you for your time, Henrik.May it be done like this above link?
Answer
Posted 9 months ago
 Henrik, If showing 11325 edges be meaningless, then is it possible to show some parts of the network like the below figure?
Answer
Posted 9 months ago
 That's great. Thank you so much, Rohit and Henrik.I found an interesting image in which the relationship between points is defined based on the coefficient of correlation between points. Can this method be evaluated for the below points based on the strength values? Subgraph[g, {"ZONGULDAK", "YUNAK", "İSKENDERUN", "BALIKESİRGÖNEN", "ÇORUM", "MUĞLA", "KONYAHAVALİMANI", "IĞDIR"}, VertexSize -> Tiny, EdgeStyle -> Arrowheads[.015]] 
Answer
Posted 9 months ago
 Dear Rohit,I meant defining a threshold like the mean strength. In other words, connecting between pairs of stations with strength more than 11.37.
Answer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments