Commit 9c004dde authored by Anthony Schrapffer's avatar Anthony Schrapffer
Browse files

Add RoutingPPViewer

parent d77c565f
#!/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()
#!/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
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)
################################################