Plotting Voronoi diagram for a country

[ghorbani (MA Ghorbani)]

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

[dauphinester (André Dauphiné)]

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

[henrikschachner]

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

The respective discussion can be found here:

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[Dynam…

[ghorbani (MA Ghorbani)]

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.

Turkey Voronoi.nb (192 KB)

[rohitnamjoshi]

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

[henrikschachner]

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]

[aliteimor (Alex Teymouri)]

Dear Henrik and Rohit,

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

Thank you

[henrikschachner]

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!
[1]: https://community.wolfram.com//c/portal/getImageAttachment?filename=rain3.png&userId=32203

[ghorbani (MA Ghorbani)]

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

I deeply appreciate your great help.

[henrikschachner]

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

[ghorbani (MA Ghorbani)]

Thanks, Henrik.

That is a perfect suggestion.

[ghorbani (MA Ghorbani)]

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,

[henrikschachner]

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

[ghorbani (MA Ghorbani)]

Dear Henrik,

That’s a nice solution. I really appreciate your time

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

[henrikschachner]

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

[ghorbani (MA Ghorbani)]

Dear Henrik,
I got this kind of output:

[henrikschachner]

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

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

Turkey_PolyArea_0.nb (22.1 KB)

[ghorbani (MA Ghorbani)]

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

data.nb (162 KB)

[henrikschachner]

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

[ghorbani (MA Ghorbani)]

Thank you for your time, Henrik.

May it be done like this above link?