From 02825e89aa60b0ef4f1481f4738b6d38b63e1fe5 Mon Sep 17 00:00:00 2001
From: Jan Polcher <jan.polcher@lmd.jussieu.fr>
Date: Sat, 15 Jun 2019 16:11:48 +0200
Subject: [PATCH] Plot of graphs improved but still not perfect.
---
Diags/plotgraphs.py | 32 +++++++++++++++++++++++++-------
run.def | 28 ----------------------------
2 files changed, 25 insertions(+), 35 deletions(-)
delete mode 100644 run.def
diff --git a/Diags/plotgraphs.py b/Diags/plotgraphs.py
index c3d3af5..b705c2f 100644
--- a/Diags/plotgraphs.py
+++ b/Diags/plotgraphs.py
@@ -4,6 +4,7 @@ import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
from netCDF4 import Dataset
+import sys
NCfile="../MEDCORDEX_test_graph.nc"
GraphicFile = "test.png"
@@ -39,6 +40,12 @@ lat_bnd = nc.variables["lat_bnd"][:,:,:]
land = nc.variables["land"][:,:]
routetogrid = nc.variables["routetogrid"][:,:,:]
routetohtu = nc.variables["routetobasin"][:,:,:]
+#
+# Coding for routetobasin :
+# nbasmax + 1 = Returnflow to the grid
+# nbasmax + 2 = Coastal flow
+# nbasmax + 3 = River flow
+#
cglon = nc.variables["CG_lon"][:,:,:]
cglat = nc.variables["CG_lat"][:,:,:]
@@ -46,7 +53,7 @@ nc.close()
nhtu, nj, ni = routetogrid.shape
landindex = LandGrid(land)
-print "Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu
+print("Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu)
#
# Get some average size of the grid boxes.
@@ -75,6 +82,8 @@ m.drawmeridians(np.arange(-80.,80.,0.1),labels=[0,0,0,1]) # draw meridians
cmap = 'RdBu'
cmap = 'Pastel1'
cmap = 'prism'
+lwidth=0.0005
+markersize=5
#
#
#
@@ -90,24 +99,33 @@ for i in range(ni) :
#
for il in range(landindex.nbland) :
j,i = landindex.land2ij(il)
- print il,"== routetogrid ==",routetogrid[:,j,i]
- print il,"== routetohtu ==",routetohtu[:,j,i]
+ xt,yt=m(lon_full[j,i],lat_full[j,i])
+ plt.text(xt,yt,str(il))
+ #
for ih in range(nhtu) :
lo = cglon[ih,j,i]
la = cglat[ih,j,i]
xl,yl=m(lo, la)
if int(routetohtu[ih,j,i]-1) < nhtu :
- plt.plot(xl,yl, marker=11)
+ plt.plot(xl,yl, ms=markersize, marker=11)
+ elif int(routetohtu[ih,j,i]-1) == nhtu+1 :
+ plt.plot(xl,yl, ms=markersize, marker="x")
+ elif int(routetohtu[ih,j,i]-1) == nhtu+2 :
+ plt.plot(xl,yl, ms=markersize, marker="o")
+ elif int(routetohtu[ih,j,i]-1) == nhtu+3 :
+ plt.plot(xl,yl, ms=markersize, marker="p")
else :
- plt.plot(xl,yl, marker="o")
+ print("Unforseen coding of outlof HTU : ",int(routetohtu[ih,j,i]-1))
+ sys.exit()
+ #
if int(routetohtu[ih,j,i]-1) < nhtu :
ile = int(routetogrid[ih,j,i]-1)
- ie,je = landindex.land2ij(ile)
+ je,ie = landindex.land2ij(ile)
ihe = int(routetohtu[ih,j,i]-1)
loe = cglon[ihe,je,ie]
lae = cglat[ihe,je,ie]
xle,yle=m(loe, lae)
- plt.plot([xl,xle],[yl,yle], color="b")
+ plt.arrow(xl, yl, xle-xl, yle-yl, color="b", width=lwidth, head_width=20*lwidth, length_includes_head=True)
#
#
#
diff --git a/run.def b/run.def
deleted file mode 100644
index 1cca8b5..0000000
--- a/run.def
+++ /dev/null
@@ -1,28 +0,0 @@
-[OverAll]
-#
-#
-EarthRadius = 6370000.
-#
-ModelGridFile = /home/polcher/WORK/Data/NewRouting/geo_em.d01.nc
-# Mallorca
-WEST_EAST = 2.3, 3.5
-SOUTH_NORTH = 39.00, 40.1
-HydroFile = /home/polcher/WORK/Data/NewRouting/routing_MED.nc
-#
-# FORTRAN interface parameters
-#
-Documentation = true
-#
-# Configuration for the graph to be generated
-#
-nbasmax = 35
-#
-# Output
-#
-GraphFile = MEDCORDEX_test_graph.nc
-#
-# Diagnostics
-# You need to provide an interval in longitude and Latitude.
-#
-DiagLon = 2.3, 3.5
-DiagLat = 39.0, 40.1
--
GitLab