From 9c004dde5ab90826019689208da6130b2c7f0263 Mon Sep 17 00:00:00 2001
From: Anthony <schrapff.ant@gmail.com>
Date: Thu, 7 May 2020 09:01:18 -0300
Subject: [PATCH] Add RoutingPPViewer
---
Diags/RoutingPPViewer/Voronoi_grid.py | 134 +++++++
Diags/RoutingPPViewer/class_Routing.py | 148 +++++++
Diags/RoutingPPViewer/interactive.py | 523 +++++++++++++++++++++++++
3 files changed, 805 insertions(+)
create mode 100644 Diags/RoutingPPViewer/Voronoi_grid.py
create mode 100644 Diags/RoutingPPViewer/class_Routing.py
create mode 100644 Diags/RoutingPPViewer/interactive.py
diff --git a/Diags/RoutingPPViewer/Voronoi_grid.py b/Diags/RoutingPPViewer/Voronoi_grid.py
new file mode 100644
index 0000000..a680363
--- /dev/null
+++ b/Diags/RoutingPPViewer/Voronoi_grid.py
@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Apr 30 13:18:21 2020
+
+@author: anthony
+"""
+
+from scipy.spatial import Voronoi, voronoi_plot_2d
+import numpy as np
+import shapely.geometry
+import shapely.ops as ops
+import shapely
+import matplotlib.pyplot as plt
+
+def voronoi_finite_polygons_2d(vor, radius=None):
+ """
+ Reconstruct infinite voronoi regions in a 2D diagram to finite
+ regions.
+
+ Parameters
+ ----------
+ vor : Voronoi
+ Input diagram
+ radius : float, optional
+ Distance to 'points at infinity'.
+
+ Returns
+ -------
+ regions : list of tuples
+ Indices of vertices in each revised Voronoi regions.
+ vertices : list of tuples
+ Coordinates for revised Voronoi vertices. Same as coordinates
+ of input vertices, with 'points at infinity' appended to the
+ end.
+
+ """
+
+ if vor.points.shape[1] != 2:
+ raise ValueError("Requires 2D input")
+
+ new_regions = []
+ new_vertices = vor.vertices.tolist()
+
+ center = vor.points.mean(axis=0)
+ if radius is None:
+ radius = vor.points.ptp().max()
+
+ # Construct a map containing all ridges for a given point
+ all_ridges = {}
+ for (p1, p2), (v1, v2) in zip(vor.ridge_points, vor.ridge_vertices):
+ all_ridges.setdefault(p1, []).append((p2, v1, v2))
+ all_ridges.setdefault(p2, []).append((p1, v1, v2))
+
+ # Reconstruct infinite regions
+ for p1, region in enumerate(vor.point_region):
+ vertices = vor.regions[region]
+
+ if all(v >= 0 for v in vertices):
+ # finite region
+ new_regions.append(vertices)
+ continue
+
+ # reconstruct a non-finite region
+ ridges = all_ridges[p1]
+ new_region = [v for v in vertices if v >= 0]
+
+ for p2, v1, v2 in ridges:
+ if v2 < 0:
+ v1, v2 = v2, v1
+ if v1 >= 0:
+ # finite ridge: already in the region
+ continue
+
+ # Compute the missing endpoint of an infinite ridge
+
+ t = vor.points[p2] - vor.points[p1] # tangent
+ t /= np.linalg.norm(t)
+ n = np.array([-t[1], t[0]]) # normal
+
+ midpoint = vor.points[[p1, p2]].mean(axis=0)
+ direction = np.sign(np.dot(midpoint - center, n)) * n
+ far_point = vor.vertices[v2] + direction * radius
+
+ new_region.append(len(new_vertices))
+ new_vertices.append(far_point.tolist())
+
+ # sort region counterclockwise
+ vs = np.asarray([new_vertices[v] for v in new_region])
+ c = vs.mean(axis=0)
+ angles = np.arctan2(vs[:,1] - c[1], vs[:,0] - c[0])
+ new_region = np.array(new_region)[np.argsort(angles)]
+
+ # finish
+ new_regions.append(new_region.tolist())
+
+ return new_regions, np.asarray(new_vertices)
+
+
+def get_voronoi_polygons(P, border):
+ vor = Voronoi(P)
+ regions, vertices = voronoi_finite_polygons_2d(vor)
+
+ polygons = []
+ for region in regions:
+ polygon = vertices[region]
+ if shapely.geometry.Polygon(polygon).intersects(border):
+ polygons.append(shapely.geometry.Polygon(polygon).intersection(border))
+ else:
+ polygons.append([])
+
+ # Equivalence points de P et polygons -> get pts de region
+ []
+
+ return polygons
+
+
+if __name__ == "__main__":
+
+ border = shapely.geometry.Polygon([[0,0], [0,1], [1,1],[1,0],[0,0]])
+ P = np.random.random((35,2))
+
+ polygons = get_voronoi_polygons(P, border)
+ x,y = border.exterior.xy
+ plt.plot(x,y)
+
+ for p in polygons:
+ p1 = p.intersection(border)
+ x,y = p1.exterior.xy
+ plt.plot(x,y)
+
+ plt.show()
+
+
diff --git a/Diags/RoutingPPViewer/class_Routing.py b/Diags/RoutingPPViewer/class_Routing.py
new file mode 100644
index 0000000..dcaee35
--- /dev/null
+++ b/Diags/RoutingPPViewer/class_Routing.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+class Routing File
+@author: anthony
+"""
+
+from netCDF4 import Dataset as NetCDFFile
+import numpy as np
+import numpy.ma as ma
+
+import matplotlib.pyplot as plt
+import cartopy
+import cartopy.crs as ccrs
+import cartopy.io.shapereader as shpreader
+from shapely.geometry import Polygon
+
+
+
+class routing:
+ def __init__(self, dfile):
+ self.nc = NetCDFFile(dfile, "r")
+ self.varlist = list(self.nc.variables.keys())
+
+ self.get_vars()
+ self.get_extent()
+ self.get_grid()
+ self.get_basin_list()
+
+ ##############################
+ # Initialization Methods
+ #
+ def get_vars(self):
+ self.var = {}
+ for v in self.varlist:
+ self.var[v] = self.nc.variables[v][:]
+ self.var["basinid"] = ma.masked_where(np.isnan(self.var["basinid"]), self.var["basinid"])
+ self.nbpt = int(np.sum(self.var["land"]))
+ self.nbasmax = self.var["basinid"].shape[0]
+ self.nbnd = self.var["lon_bnd"].shape[0]
+
+ self.conv_land2ij = [[np.where(self.var["nbpt_glo"] == i)[0][0], \
+ np.where(self.var["nbpt_glo"] == i)[1][0]] \
+ for i in range(1, self.nbpt+1)]
+ self.var_land = {}
+ for v in self.varlist:
+ if self.var[v].ndim == 2:
+ self.var_land[v] = [self.var[v][j,i] for j,i in self.conv_land2ij]
+ elif v == "lon_bnd" or v == "lat_bnd":
+ self.var_land[v] = [self.var[v][:,j,i] \
+ for j,i in self.conv_land2ij]
+ elif self.var[v].ndim == 3:
+ self.var_land[v] = [self.var[v][:int(self.var["routenbintobas"][j,i]),j,i] \
+ for j,i in self.conv_land2ij]
+
+ def get_grid(self):
+ self.GRID = [self.get_polygon(i) for i in range(self.nbpt)]
+ def get_polygon(self, i):
+ pts = [(self.var_land["lon_bnd"][i][l], self.var_land["lat_bnd"][i][l]) for l in range(self.nbnd)]
+ return Polygon(pts)
+
+ def get_extent(self, bid = None):
+ self.resh = np.nanmax(np.diff(self.var["lon"][:,:], axis = 1))
+ self.resv = np.nanmax(np.diff(self.var["lat"][:,:], axis = 0))
+
+ if bid is not None:
+ mask = self.basin_mask(bid)
+ mask = np.sum(mask, axis = 0)
+ A = np.where(mask>0)
+ lonnd = [self.var["lon_bnd"][:,A[0][i], A[1][i]] for i in range(len(A[0]))]
+ latnd = [self.var["lat_bnd"][:, A[0][i], A[1][i]] for i in range(len(A[0]))]
+
+ ext = [np.min(lonnd)-self.resh, np.max(lonnd)+self.resh, \
+ np.min(latnd)-self.resv, np.max(latnd)+self.resh]
+ return ext
+ else:
+ self.ext = [np.nanmin(self.var["lon_bnd"])-self.resh, np.nanmax(self.var["lon_bnd"])+self.resh, \
+ np.nanmin(self.var["lat_bnd"])-self.resv, np.nanmax(self.var["lat_bnd"])+self.resv]
+ ##############################
+ # Basic basin operation
+ #
+ def get_basin_list(self):
+ basin_list = ma.masked_where(np.isnan(self.var["basinid"]),self.var["basinid"])
+ basin_list = list(ma.unique(basin_list[basin_list.mask == False]))
+ area = [np.sum(self.basin_area(bid)) for bid in basin_list]
+ self.basin_list = [[int(basin_list[i]), int(area[i])] for i in np.flip(np.argsort(area))]
+
+ def basin_area(self,bid):
+ mask = self.basin_mask(bid)
+ area = np.nansum(mask*self.var["basin_area"], axis = 0)
+ return area
+
+ def basin_mask(self,bid):
+ A = ma.where(self.var["basinid"] == float(bid))
+ pts = self.where2list(A)
+ #
+ mask = np.zeros(self.var["basinid"].shape)
+ mask[A] = 1
+ return mask
+
+ def basin_listpoints(self, bid):
+ mask = np.sum(self.basin_mask(bid), axis = 0)
+ lstpts = [i for i in range(self.nbpt) if mask[r1.conv_land2ij[i][0],r1.conv_land2ij[i][1]]> 0]
+ return lstpts
+
+ def main_fetch(self, bid, p = 90):
+ mask = self.basin_mask(bid)
+ pts = self.where2list(np.where(mask==1))
+ fetch = [self.var["fetch"][a,b,c] for a,b,c in pts]
+ #
+ maxfetch = np.max(np.array(fetch))
+ perc = np.percentile(fetch, p)
+ A = np.where(fetch>perc)[0]
+ pts_out = [pts[i] for i in A]
+ return pts_out
+
+ def outflow(self, bid):
+ mask = self.basin_mask(bid)
+ pts = self.where2list(np.where(mask==1))
+ out1 = ma.filled(self.var["routetobasin"], -99)
+ outbas1 = np.array([out1[a,b,c] for a,b,c in pts])
+ pts_out = [pts[i]for i in np.where(outbas1>self.nbasmax)[0]]
+ outflow = [[int(self.var["routetobasin"][a,b,c]), int(self.var["fetch"][a,b,c]),[a,b,c]] for a,b,c in pts_out]
+
+ return outflow
+
+ def inflow_outflow_htus(self, htus):
+ x = []
+ y = []
+ for a,b,c in htus:
+ k = int(self.var["routetobasin"][a,b,c]-1)
+ if k<self.nbasmax:
+ jj,ii = self.conv_land2ij[int(self.var["routetogrid"][a,b,c]-1)]
+ x.append((self.var["CG_lon"][a,b,c], self.var["CG_lon"][k,jj,ii]))
+ y.append((self.var["CG_lat"][a,b,c], self.var["CG_lat"][k,jj,ii]))
+ else:
+ x.append([self.var["CG_lon"][a,b,c]])
+ y.append([self.var["CG_lat"][a,b,c]])
+
+ return x,y
+
+ ##############################
+ # Basic Methods
+ #
+ def where2list(self, A):
+ pts = [[A[0][i],A[1][i], A[2][i]] for i in range(len(A[0]))]
+ return pts
+
diff --git a/Diags/RoutingPPViewer/interactive.py b/Diags/RoutingPPViewer/interactive.py
new file mode 100644
index 0000000..a6a4f07
--- /dev/null
+++ b/Diags/RoutingPPViewer/interactive.py
@@ -0,0 +1,523 @@
+from tkinter import *
+from tkinter import font
+from tkinter import filedialog
+import sys
+import matplotlib
+matplotlib.use("TkAgg")
+from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+import cartopy.io.shapereader as shpreader
+import cartopy.io.img_tiles as cimgt
+import numpy as np
+import matplotlib as mpl
+from matplotlib import cm
+
+# Add the location where the programs are located if necessary
+#sys.path.append("/home/anthony/Documents/LOCAL/RoutingPP/ROUTING_Analysis_tools/Programs")
+
+from class_Routing import routing
+from Voronoi_grid import get_voronoi_polygons
+import numpy.ma as ma
+
+################################################
+#
+# Class Interface
+# Main class of the program
+
+class Interface(Frame):
+ def __init__(self, parent, **kwargs):
+ Frame.__init__(self, parent, width=768, height=576)
+
+ self.container = Frame(self)
+ self.container.pack(side="top", fill="both", expand = True)
+ self.container.grid_rowconfigure(0, weight=1)
+ self.container.grid_columnconfigure(0, weight=1)
+ self.parent = parent
+ self.loaded = False
+ self.grid = False
+ self.main = False
+ self.basin_id = 0
+ self.nbpt = 0
+ self.HTUview = "global" # or local
+ self.HTUneighbours = False
+ self.HTUarrow = False
+ self.initUI()
+
+
+ def initUI(self):
+ self.parent.title("RoutingPP Viewer")
+ self.pack(fill=BOTH, expand=1)
+
+ menubar = Menu(self.parent)
+ self.parent.config(menu=menubar)
+
+ fileMenu = Menu(menubar)
+ fileMenu.add_command(label="Open", command=self.onOpen)
+ fileMenu.add_command(label="Exit", command=self.parent.quit)
+
+ menubar.add_cascade(label="File", menu=fileMenu)
+
+ self.frames = {}
+ frame = BASINS(self.container, self)
+ frame.grid(row=0, column=0, sticky="nsew")
+ self.frames[BASINS] = frame
+ self.show_frame(BASINS)
+
+ def show_frame(self, cont):
+ frame = self.frames[cont]
+ frame.grid(row=0,column=0,sticky="nsew")
+ frame.tkraise()
+
+ def update_BASINS(self):
+ frame = BASINS(self.container, self, self.title)
+ frame.grid(row=0, column=0, sticky="nsew")
+ self.frames[BASINS] = frame
+ self.show_frame(BASINS)
+
+ def update_HTUS(self):
+ frame = HTUS(self.container, self, self.title)
+ frame.grid(row=0, column=0, sticky="nsew")
+ self.frames[HTUS] = frame
+ self.show_frame(HTUS)
+
+ def get_nbpt(self, nbpt):
+ self.nbpt = nbpt
+ self.update_HTUS()
+
+ def grid_button(self):
+ self.grid = not self.grid
+ self.update_BASINS()
+
+ def basin_plot(self):
+ self.update_BASINS()
+
+ def main_river(self):
+ self.main = not self.main
+ self.update_BASINS()
+
+ def local_global(self):
+ if self.HTUview == "global" :
+ self.HTUview = "local"
+ elif self.HTUview == "local" :
+ self.HTUview = "global"
+ self.update_HTUS()
+
+ def arrows(self):
+ self.HTUarrow = not self.HTUarrow
+ self.update_HTUS()
+
+ def neighbours(self):
+ self.HTUneighbours = not self.HTUneighbours
+ self.update_HTUS()
+
+ def onOpen(self):
+
+ ftypes = [('NetCDF Routing', '*.nc'), ('All files', '*')]
+ dlg = filedialog.Open(self, filetypes = ftypes)
+ fl = dlg.show()
+
+ # to test the program
+ #fl = "/home/anthony/Documents/LOCAL/RoutingPP/ROUTING_Analysis_tools/Originals/LC_Spain_test_graph_corr.nc"
+ #fl = "/home/anthony/Documents/LOCAL/RoutingPP/ROUTING_Analysis_tools/Originals/EuroSW_Mallorca.nc"
+
+ self.frames = {}
+ if fl != '':
+ self.rr = routing(fl)
+ self.loaded = True
+ self.title = fl.split("/")[-1]
+ frame = BASINS(self.container, self, self.title)
+ frame.grid(row=0, column=0, sticky="nsew")
+ self.frames[BASINS] = frame
+ self.show_frame(BASINS)
+
+################################################
+#
+# Class BASIN
+# Class of the page showing the basins
+
+class BASINS(Frame):
+ def __init__(self, parent, controller, title = "River routing, select a file"):
+ Frame.__init__(self, parent)
+ label = Label(self, text=title)
+ label.pack(pady=10,padx=10)
+ self.controller = controller
+
+ if controller.loaded:
+
+ b1 = Label(self, text="Actual mode : Basin view")
+ b1.pack(pady=10,padx=10)
+ #b1.pack(side = TOP)
+
+ b2 = Button(self, text="Switch to HTUs view", command = controller.update_HTUS)
+ b2.pack(side = TOP)
+
+ self.main_fig(controller.rr.ext)
+
+ button1 = Button(self, text="Grid", command = controller.grid_button)
+ button1.pack(side = TOP)
+
+
+ readOnlyText = Text(self)
+ label1 = Label(self, text =" Select basin : ")
+ label1.pack(side = TOP)
+
+ OPTIONS = ["None"] + [b[0] for b in controller.rr.basin_list[0:30]]
+ self.variable = StringVar(parent)
+ if controller.basin_id == 0:
+ self.variable.set("None")
+ else:
+ self.variable.set(controller.basin_id)
+ w = OptionMenu(self, self.variable, *OPTIONS)
+ w.pack(side = TOP)
+
+ button2 = Button(self, text="Show basin", command = self.basin_id)
+ button2.pack(side = TOP)
+
+ if controller.basin_id != 0:
+ button3 = Button(self, text="Main river", command = controller.main_river)
+ button3.pack(side = TOP)
+
+ else:
+ self.main_fig()
+
+ def basin_id(self):
+ if self.variable.get() == "None":
+ self.controller.basin_id = 0
+ else:
+ self.controller.basin_id = self.variable.get()
+ self.controller.basin_plot()
+
+ def main_fig(self, extent = None):
+
+ df = draw_Figure(self, extent, self.controller.grid, self.controller.basin_id, self.controller.main)
+
+ canvas = FigureCanvasTkAgg(df.fig, master=self)
+ canvas.draw()
+ canvas.get_tk_widget().pack(side=LEFT, fill=BOTH, expand=True)
+ canvas._tkcanvas.pack(side=LEFT, fill=BOTH, expand=True)
+
+################################################
+#
+# Class HTUs
+# Class of the page showing the basins
+
+class HTUS(Frame):
+
+ def __init__(self, parent, controller, title = ""):
+ Frame.__init__(self, parent)
+ label = Label(self, text=title)
+ label.pack(pady=10,padx=10)
+ self.controller = controller
+
+ b1 = Label(self, text="Actual mode : HTUs View", command = None)
+ b1.pack(pady=10,padx=10)
+
+ b2 = Button(self, text="Switch to Basin view", command = controller.update_BASINS)
+ b2.pack(side = TOP)
+
+ if controller.HTUview == "global": # or local
+ if controller.nbpt == 0:
+ self.main_fig(controller.rr.ext)
+ if controller.nbpt != 0:
+ self.main_fig(controller.rr.ext, True)
+ elif controller.HTUview == "local":
+ self.voronoi()
+
+
+ # Pour entrer le nbpt !
+ b1 = Label(self, text="Enter the grid point number")
+ b1.pack(pady=10,padx=10)
+ b1 = Label(self, text="(min: 1; max {0})".format(controller.rr.nbpt))
+ b1.pack(pady=10,padx=10)
+ self.e = Entry(self)
+ self.e.pack(pady=10,padx=10)
+
+ button1 = Button(self, text="Show selected grid point", command = self.get_nbpt)
+ button1.pack(side = TOP)
+
+ if controller.HTUview == "local":
+ text_button4 = "Switch to global map"
+
+ if controller.HTUneighbours:
+ text_button2 = "Hide neighbours"
+ else:
+ text_button2 = "Show neighbours"
+
+ button2 = Button(self, text=text_button2, command = self.neighbours)
+ button2.pack(side = TOP)
+
+ if controller.HTUarrow:
+ text_button3 = "Hide arrows"
+ else:
+ text_button3 = "Show arrows"
+
+ button3 = Button(self, text=text_button3, command = self.arrows)
+ button3.pack(side = TOP)
+ else:
+ text_button4 = "Switch to local map"
+ button4 = Button(self, text=text_button4, command = self.local_global_map)
+ button4.pack(side = TOP)
+
+ def get_nbpt(self):
+ nbpt = self.e.get()
+ if nbpt.isnumeric():
+ nbpt = int(nbpt)
+ if (nbpt>0) and (nbpt<=self.controller.rr.nbpt):
+ self.controller.get_nbpt(nbpt)
+ else:
+ print("error")
+ else:
+ print("error")
+ # error
+ return
+
+ def local_global_map(self):
+ self.controller.local_global()
+
+ def neighbours(self):
+ self.controller.neighbours()
+
+ def arrows(self):
+ self.controller.arrows()
+
+ def main_fig(self, extent = None, grid = False):
+ df = draw_Figure(self, extent, grid, pt = self.controller.nbpt)
+
+ canvas = FigureCanvasTkAgg(df.fig, master=self)
+ canvas.draw()
+ canvas.get_tk_widget().pack(side=LEFT, fill=BOTH, expand=True)
+ canvas._tkcanvas.pack(side=LEFT, fill=BOTH, expand=True)
+
+ def voronoi(self):
+ df = draw_Voronoi(self, nb = self.controller.nbpt, neighbour = self.controller.HTUneighbours, arrow = self.controller.HTUarrow)
+
+ canvas = FigureCanvasTkAgg(df.fig, master=self)
+ canvas.draw()
+ canvas.get_tk_widget().pack(side=LEFT, fill=BOTH, expand=True)
+ canvas._tkcanvas.pack(side=LEFT, fill=BOTH, expand=True)
+
+################################################
+#
+# Class draw_figure
+# Class of the map figures
+
+
+class draw_Figure:
+ """
+ Global map and local map of the basins
+ """
+ def __init__(self, page, extent, grid, basin_id = 0 , main=False, pt = 0):
+ self.page = page
+ self.extent = extent
+ self.grid = grid
+ self.basin_id = basin_id
+ self.main = main
+ self.pt = pt
+
+ self.init_figure()
+ if self.extent is not None:
+ self.local_map()
+ self.draw_grid()
+ self.draw_basin()
+ self.main_river()
+ self.draw_nbpt_loc()
+ else:
+ self.global_map()
+
+ #####
+
+ def init_figure(self):
+ self.fig = plt.figure(figsize=(8,4), dpi=100)
+ self.ax = self.fig.add_axes([0.01, 0.01, 0.98, 0.98], projection=ccrs.PlateCarree())
+
+ def draw_grid(self):
+ if self.grid:
+ for l in self.page.controller.rr.GRID:
+ x, y = l.exterior.xy
+ plt.plot(x,y, color = "k")
+
+ def draw_basin(self):
+ if self.basin_id != 0 :
+ area = self.page.controller.rr.basin_area(int(self.basin_id))
+ area = [area[j,i]/self.page.controller.rr.var["area"][j,i]*100 for j,i in self.page.controller.rr.conv_land2ij]
+ cmap = cm.get_cmap('viridis', 12)
+ norm = mpl.colors.Normalize(vmin=0.,vmax=np.max(area))
+ for l, ar in zip(self.page.controller.rr.GRID, area):
+ if ar>0:
+ x, y = l.exterior.xy
+ plt.fill(x,y, color = cmap(norm(ar)))
+ plt.title("basin_id : {0}".format(int(self.basin_id)))
+
+ def draw_nbpt_loc(self):
+ if self.pt != 0:
+ l = self.page.controller.rr.GRID[self.pt-1]
+ x, y = l.exterior.xy
+ plt.fill(x,y, color = "r")
+
+ def global_map(self):
+ self.ax.set_global()
+ self.ax.stock_img()
+ self.ax.coastlines()
+ lw = 0.5
+ river_50m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m', facecolor='none', edgecolor = "#011f4b", lw = lw)
+ self.ax.add_feature(river_50m)
+
+ def local_map(self):
+ self.ax.set_extent(self.extent)
+ if (self.basin_id != 0) or (self.pt != 0):
+ fc_coast = "none"
+ fc_ocean = "none"
+ else:
+ fc_coast ='#77ab59'
+ fc_ocean ='#99ccff'
+ coast = cfeature.NaturalEarthFeature('physical', 'coastline', '50m', facecolor=fc_coast, edgecolor = "k", lw = 1)
+ self.ax.add_feature(coast)
+ ocean = cfeature.NaturalEarthFeature('physical', 'ocean', '50m', facecolor=fc_ocean, edgecolor = "none")
+ self.ax.add_feature(ocean)
+ lw = 1
+ river_50m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m', facecolor='none', edgecolor = "#011f4b", lw = lw)
+ self.ax.add_feature(river_50m)
+
+
+ def main_river(self):
+ if self.main:
+ pts_out = self.page.controller.rr.main_fetch(int(self.basin_id), p = 90)
+ outflow = self.page.controller.rr.outflow(int(self.basin_id))
+ x,y = self.page.controller.rr.inflow_outflow_htus(pts_out)
+
+ i=0
+ for xx, yy in zip(x,y):
+ if len(xx)>1:
+ plt.plot(xx,yy,ls = "-", color = "k")
+
+ for xx, yy in zip(x,y):
+ plt.plot(xx[0], yy[0], "bo", ms = 5)
+
+ x,y = self.page.controller.rr.inflow_outflow_htus([outflow[i][2] for i in range(len(outflow))])
+ for xx, yy in zip(x,y):
+ plt.plot(xx[0], yy[0], "ro", ms = 10)
+
+ plt.title("Basin {0} main htus".format(self.basin_id))
+
+
+################################################
+#
+# Class draw_figure
+# Class of the Voronoi diagram figures
+
+class draw_Voronoi:
+ """
+ Global map and local map of the basins
+ """
+ def __init__(self, page, nb, neighbour, arrow):
+
+ if neighbour:
+ j,i = page.controller.rr.conv_land2ij[nb-1]
+ array = ma.array(page.controller.rr.var["nbpt_glo"][j-1:j+2, i-1:i+2]).astype(np.int32)
+ gridpts = list(ma.unique(array[(array.mask == False)]))
+ LNB = [pt for pt in gridpts if pt>0]
+ else:
+ LNB = [nb]
+
+ self.page = page
+ self.start_map()
+ self.get_limits(LNB)
+
+ for pt in LNB:
+ main = False
+ if pt == nb: main = True
+ self.plot_voronoi(pt, arrow, main)
+
+ x, y = page.controller.rr.GRID[nb-1].exterior.xy
+ plt.plot(x,y, color = "r")
+
+ cbaxes = self.fig.add_axes([0.82, 0.1, 0.03, 0.75])
+ cb1 = mpl.colorbar.ColorbarBase(ax = cbaxes, cmap=self.viridis,
+ norm=self.norm,
+ orientation='vertical', alpha = 0.6)
+ cb1.set_label('\% of total area of the grid box')
+ self.fig.subplots_adjust(left = 0.05, right = 0.8)
+
+
+ def start_map(self):
+ self.viridis = cm.get_cmap('Reds', 12)
+ self.norm = mpl.colors.Normalize(vmin=0.,vmax=100)
+
+ self.fig = plt.figure(figsize= (7,7))
+ self.ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())
+ geo_reg_shp = shpreader.natural_earth(resolution='10m', category='physical',\
+ name='coastline')
+ geo_reg = shpreader.Reader(geo_reg_shp)
+ self.ax.add_geometries(geo_reg.geometries(), ccrs.PlateCarree(), \
+ edgecolor="k", facecolor='none')
+
+ def get_limits(self, LNB):
+ ext = self.page.controller.rr.resh/2.
+ lat = []
+ lon = []
+ for pt in LNB:
+ j,i = self.page.controller.rr.conv_land2ij[pt-1]
+ lat = lat + list(self.page.controller.rr.var["lat_bnd"][:,j,i])
+ lon = lon + list(self.page.controller.rr.var["lon_bnd"][:,j,i])
+
+ limits = [np.min(lon)-ext, np.max(lon)+ext, np.min(lat)-ext, np.max(lat)+ext]
+ self.ax.set_extent(limits)
+
+ def plot_voronoi(self, nb, arrow = True, main = False):
+ rr = self.page.controller.rr
+ j,i = rr.conv_land2ij[nb-1]
+ nbas = int(rr.var["routenbintobas"][j,i])
+ outgrid = [g for g in rr.var["routetogrid"][:nbas,j,i]]
+ outbas = [b for b in rr.var["routetobasin"][:nbas,j,i]]
+
+ area = np.array([b for b in rr.var["basin_area"][:nbas,j,i]])
+ area =area/rr.var["area"][j,i]*100
+
+ htus = [[lon, lat] for lon, lat in zip(rr.var["CG_lon"][:nbas,j,i], rr.var["CG_lat"][:nbas,j,i])]
+ htus = np.array(htus)
+
+ X = []
+ Y = []
+ for l, htu in enumerate(htus):
+ jj,ii = rr.conv_land2ij[int(rr.var["routetogrid"][l,j,i]-1)]
+ k = int(rr.var["routetobasin"][l,j,i]-1)
+ if (k<rr.nbasmax):
+ x = [htu[0], rr.var["CG_lon"][k,jj,ii]]
+ y = [htu[1], rr.var["CG_lat"][k,jj,ii]]
+ X.append(x)
+ Y.append(y)
+
+ if arrow:
+ for x, y in zip(X, Y):
+ plt.arrow(x[0], y[0], x[1]-x[0], y[1]-y[0], "k", length_includes_head = True, head_length = 0.01, head_width = 0.01)
+
+ A, B = htus[:,0], htus[:,1]
+ plt.scatter(A, B, marker='.')
+
+ polygons = get_voronoi_polygons(htus, rr.GRID[nb-1])
+ for i,p in enumerate(polygons):
+ if p != []:
+ x,y = p.exterior.xy
+ plt.fill(x,y, color = self.viridis(self.norm(area[i])), alpha = 0.6)
+ plt.plot(x,y,color = "k")
+
+
+ x, y = rr.GRID[nb-1].exterior.xy
+ plt.plot(x,y, color = "k")
+
+ i = 0
+ for lon, lat in htus:
+ plt.plot(lon, lat, "o",alpha = 0.7, color = "b")
+ if outbas[i]>rr.nbasmax:
+ plt.plot(lon, lat, "o",alpha = 0.7, ms = 5, color = "r")
+ i+=1
+
+def main():
+ root = Tk()
+ interface = Interface(root)
+ root.mainloop()
+
+if __name__ == '__main__':
+ main()
--
GitLab