HydroGrid.py 6.44 KB
Newer Older
1 2 3 4 5 6 7 8
#
# Reads ans manages the model grid
#
import numpy as np
from netCDF4 import Dataset
# Doc : https://spacetelescope.github.io/spherical_geometry/api/spherical_geometry.polygon.SphericalPolygon.html
from spherical_geometry import polygon
import RPPtools as RPP
9
import os 
10 11 12 13
#
import configparser
config=configparser.ConfigParser()
config.read("run.def")
14
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
15 16 17 18 19 20 21
#
import getargs
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
#
def getbox(ncdf, corners) :
22 23 24 25 26 27
    # Add a few points to tbe box to make sure to cover everything
    halo_pts=2
    ii=[np.argmin(abs(ncdf.variables["nav_lon"][0,:]-np.min(corners[0][:]))),
        np.argmin(abs(ncdf.variables["nav_lon"][0,:]-np.max(corners[0][:])))]
    jj=[np.argmin(abs(ncdf.variables["nav_lat"][:,0]-np.min(corners[1][:]))),
        np.argmin(abs(ncdf.variables["nav_lat"][:,0]-np.max(corners[1][:])))]
28

29
    return min(ii)-halo_pts, max(ii)+halo_pts, min(jj)-halo_pts, max(jj)+halo_pts
30 31 32
#
# Get the corners for a regular lat/lon grid
#
33
def corners(lon, lat, hdlon, hdlat) :
34
    jjm,iim = lon.shape
35 36
    #
    index = [[j,i] for i in range(iim) for j in range(jjm)]
37 38 39 40 41
    #
    cornersll = [RPP.boxit(np.array([lon[j,i], lat[j,i]]), hdlon, hdlat, 2) for j,i in index]
    Llats = [[0.0000001 if p[1] == 0 else np.sign(p[1])*89.99 if np.abs(p[1]) == 90 else p[1] for p in boxll] for boxll in cornersll]

    Llons = [[0.0000001 if p[0] == 0 else np.sign(p[0])*89.99 if np.abs(p[0]) == 90 else p[0] for p in boxll] for boxll in cornersll]
42 43
    centersll = [[lon[j,i], lat[j,i]] for j,i in index]
    cornerspoly = [ polygon.SphericalPolygon.from_lonlat(lons, lats, center=cent) for lons, lats, cent in zip(Llons, Llats, centersll)]
Anthony's avatar
Anthony committed
44
    radiusll = [RPP.maxradius(np.array(cent), np.array(lons), np.array(lats)) for cent, lons, lats in zip(centersll,Llons,Llats)]
45
    #
46
    return cornersll, cornerspoly, centersll, radiusll, index
47
#
48
def gather(x, index, default = 0) :
49 50
    y=[]
    for ia in index :
51 52 53 54
        if (ia[:,0] == [-1, -1]).all():
            y.append([default])
        else:
            y.append(list(x[ia[0,i],ia[1,i]] for i in range(ia.shape[1]) ))
55 56 57 58 59 60 61 62 63 64 65
    return y
#
def getattrcontaining(nc, varname, substr) :
    att=[]
    for s in nc.variables[varname].ncattrs() :
        if s.lower().find(substr) >= 0 :
            att.append(nc.variables[varname].getncattr(s))

    return att
#
class HydroGrid :
66
    def __init__(self, lolacorners, wfile) :
67
        #
68 69 70
        self.source=config.get("OverAll", "HydroFile")
        INFO("Opening in HydroGrid : "+self.source)
        self.ncfile=Dataset(self.source,'r')
71
        istr, iend, jstr, jend = getbox(self.ncfile, lolacorners)
72 73 74 75 76 77 78 79
        hdlon = np.mean(np.abs(np.diff(self.ncfile.variables["nav_lon"][0,:])))
        hdlat = np.mean(np.abs(np.diff(self.ncfile.variables["nav_lat"][:,0])))
        # In case the box is on the edge of the routing.nc file
        if istr<0:
            istr = 0  
        if iend >=self.ncfile.variables["nav_lon"][0,:].shape[0]: 
            iend = self.ncfile.variables["nav_lon"][0,:].shape[0]
        self.box=[istr, iend, jstr, jend]   
80 81 82 83 84
        self.lon=np.copy(self.ncfile.variables["nav_lon"][jstr:jend,istr:iend])
        self.lat=np.copy(self.ncfile.variables["nav_lat"][jstr:jend,istr:iend])
        self.jjm,self.iim = self.lon.shape
        DEBUG("# Range Lon :"+str(np.min(self.lon))+" -- "+str(np.max(self.lon)))
        DEBUG("# Range Lat :"+str(np.min(self.lat))+" -- "+str(np.max(self.lat)))
85
        #
86
        if wfile is None or not os.path.exists(wfile): 
87
            self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat, hdlon, hdlat)
88 89

    def select(self, c, r) :
Anthony's avatar
Anthony committed
90
        indices = [i for i in range(len(self.centers)) if (RPP.loladist(np.array(c),np.array(self.centers[i])) <= r+self.radius[i])]
91
        return indices
92 93 94 95

class HydroData :
    def __init__(self, nf, box, index) :
        istr, iend, jstr, jend = box[:]
96 97 98
        #
        # Flow direction
        #
99
        self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index, 97)
100 101
        self.tripdesc=nf.variables["trip"].long_name
        #
102 103
        # ID of basin
        #
104
        self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index, 999)
105 106 107 108 109 110
        self.basinsdesc=nf.variables["basins"].long_name
        att = getattrcontaining(nf, "basins", "max")
        if len(att) > 0 :
            self.basinsmax=att[0]
        else :
            INFO("We need to scan full file to find maximum number of basins")
111 112
            # This variable seems not to be used further
            self.basinsmax = part.domainmax(np.ma.max(ma.masked_where(self.basins < 1.e10, self.basins)))
113
        #
114
        # Distance to ocean
115
        #
116
        self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index, 0)
117 118
        self.distodesc=nf.variables["disto"].long_name
        #
119 120
        # Flow accumulation
        #
121
        self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index, 0)
122
        self.facdesc=nf.variables["fac"].long_name
123
        #
124 125 126
        # If Orography is present in the file then we can compute the topographic index.
        # Else the topographic index needs to be present in the file.
        #
127
        if "orog" in nf.variables.keys():
128
            self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index, 0)
129 130 131
            self.orogdesc=nf.variables["orog"].long_name
        else:
            self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
132 133 134 135 136 137 138 139 140
            self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, 10)
            self.topoinddesc=nf.variables["topoind"].long_name
            att = getattrcontaining(nf, "topoind", "min")
            if len(att) > 0 :
                self.topoindmin=att[0]
            else :
                INFO("We need to scan full file to find minimum topoind over domain")
                self.topoindmin=np.min(np.where(nf.variables["topoind"][:,:] <  1.e15))
        #
141 142
        #
        if "floodplains" in nf.variables.keys():
143
            self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index, 0)
144 145 146
            self.floodplainsdesc=nf.variables["floodplains"].long_name
        else:
            self.floodplains = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)