Message Boards Message Boards

GROUPS:

Plotting Voronoi diagram for a country

Posted 7 months ago
1992 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]}]

enter image description here

25 Replies

You have a similar question by André Dauphiné and a response by Henrik Schancher A. Dauphiné

@André Dauphiné: Thanks for the reference and greetings!

The respective discussion can be found here:

https://community.wolfram.com/groups/-/m/t/410510

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.

Attachment

Posted 7 months 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]]

enter image description here

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,

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

Dear Henrik,

That's a nice solution. I really appreciate your time :)

Is it possible to add station names to the output table?

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"}]

Dear Henrik, I got this kind of output: enter image description here

Ok, then you did not evaluate eventually a line like stations = ...

Anyway - please find attached my full notebook. Regards -- Henrik

Attachments:

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]

enter image description here

Posted 7 months ago

Dear Henrik and Rohit,

Awesome solution for a difficult question. It’s so helpful.

Thank you :)

... 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]

enter image description here

Thanks, Henrik.

That is a perfect suggestion.

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-data

For 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.

Regards

Mohammad

Attachments:

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"<s. If so, you can try:

@Moderation Team: I need Blockquote here in order to get the code below formatted correctly ...

stats = MapThread[Association["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

Thank you for your time, Henrik.

May it be done like this above link?

Posted 4 months ago

Henrik, If showing 11325 edges be meaningless, then is it possible to show some parts of the network like the below figure? enter image description here

Posted 4 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]]

enter image description here

Posted 4 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. enter image description here

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]]
Posted 4 months ago

Alex,

How is "linear correlation coefficient" defined based on strength values?

Posted 4 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. enter image description here

Posted 4 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]

enter image description here

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]

enter image description here

@Alex Teymouri : Thank you for your kind words!

Rohit and Henrik, You did a perfect and nice solution.

I deeply appreciate your great help.

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