plotgraphs.py 3.84 KB
Newer Older
1 2 3 4 5 6
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import numpy as np
from netCDF4 import Dataset
7
import sys
8

9
NCfile="../SPortugal_test_graph_n4.nc"
10
GraphicFile = "test.png"
11
title="From file : "+NCfile
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
#
#####################################################################
#
class LandGrid :
    def __init__(self, land) :
        nj,ni = land.shape
        self.indP=[]
        self.nbland = 0
        for i in range(ni) :
            for j in range(nj) :
                if (land[j,i] > 0 ) :
                    self.nbland = self.nbland+1
                    self.indP.append([j,i])
        return
                    
    def land2ij(self, il) :
        j,i = self.indP[il]
        return j,i
#
#####################################################################
#
#
nc = Dataset(NCfile, 'r')

lon_full = nc.variables["lon"][:,:]
lat_full = nc.variables["lat"][:,:]
lon_bnd = nc.variables["lon_bnd"][:,:,:]
lat_bnd = nc.variables["lat_bnd"][:,:,:]
land = nc.variables["land"][:,:]
routetogrid = nc.variables["routetogrid"][:,:,:]
routetohtu = nc.variables["routetobasin"][:,:,:]
43 44 45 46 47 48
#
# Coding for routetobasin :
# nbasmax + 1 = Returnflow to the grid
# nbasmax + 2 = Coastal flow
# nbasmax + 3 = River flow
#
49 50 51 52 53 54 55
cglon = nc.variables["CG_lon"][:,:,:]
cglat = nc.variables["CG_lat"][:,:,:]

nc.close()

nhtu, nj, ni = routetogrid.shape
landindex = LandGrid(land)
56
print("Nb land points = ", landindex.nbland, 'Nb HTU : ', nhtu)
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84

#
# Get some average size of the grid boxes.
#
lon_res = np.nanmean(np.abs(np.gradient(lon_bnd, axis=0)))
lat_res = np.nanmean(np.abs(np.gradient(lat_bnd, axis=0)))
#
#
#
fig = plt.figure(1)
ax = fig.add_subplot(111)
ax.set_title(title)
ax.set_xlabel('Longitude',labelpad=20)
ax.set_ylabel('Latitude',labelpad=20)
ax.axis([0, 10, 0, 10])
#
m = Basemap(projection='cyl', resolution='i', \
            llcrnrlon=np.nanmin(lon_bnd)-lon_res, llcrnrlat=np.nanmin(lat_bnd)-lat_res, \
            urcrnrlon=np.nanmax(lon_bnd)+lon_res, urcrnrlat=np.nanmax(lat_bnd)+lat_res)

m.drawcoastlines(linewidth=0.5)
m.drawcountries(linewidth=0.5)
m.drawparallels(np.arange(-180.,180.,0.1),labels=[1,0,0,0]) # draw parallels
m.drawmeridians(np.arange(-80.,80.,0.1),labels=[0,0,0,1]) # draw meridians
#
cmap = 'RdBu'
cmap = 'Pastel1'
cmap = 'prism'
85 86
lwidth=0.0005
markersize=5
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
#
#
#
for i in range(ni) :
    for j in range(nj) :
        if land[j,i] > 0 :
            lo = np.append(lon_bnd[:,j,i],lon_bnd[0,j,i])
            la = np.append(lat_bnd[:,j,i],lat_bnd[0,j,i])
            xl,yl=m(lo, la)
            plt.plot(xl,yl, color="m")
#
# Draw the graph
#
for il in range(landindex.nbland) :
    j,i = landindex.land2ij(il)
102 103
    xt,yt=m(lon_full[j,i],lat_full[j,i])
    plt.text(xt,yt,str(il))
104 105
    act_nhtu=np.count_nonzero(routetohtu[:,j,i] > 0)
    for ih in range(act_nhtu) :
106 107 108 109
        lo = cglon[ih,j,i]
        la = cglat[ih,j,i]
        xl,yl=m(lo, la)
        if int(routetohtu[ih,j,i]-1) < nhtu :
110 111 112 113 114 115 116
            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")
117
        else :
118 119
            print("Unforseen coding of outlof HTU : ",int(routetohtu[ih,j,i]-1))
            sys.exit()
120
        
121 122
        if int(routetohtu[ih,j,i]-1) < nhtu :
            ile =  int(routetogrid[ih,j,i]-1)
123
            je,ie = landindex.land2ij(ile)
124 125 126 127
            ihe = int(routetohtu[ih,j,i]-1)
            loe = cglon[ihe,je,ie]
            lae = cglat[ihe,je,ie]
            xle,yle=m(loe, lae)
128
            plt.arrow(xl, yl, xle-xl, yle-yl, color="b", width=lwidth, head_width=20*lwidth, length_includes_head=True)
129 130 131 132 133
#
#
#
plt.savefig(GraphicFile, dpi=450)
plt.close()