HydroGrid.py 5.78 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
#
# 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
#
import configparser
config=configparser.ConfigParser()
config.read("run.def")
EarthRadius=config.getfloat("OverAll", "EarthRadius")
#
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) :
    halo=0.1
    ii=[np.argmin(abs(ncdf.variables["nav_lon"][0,:]-(np.min(corners[0][:])-halo))),
        np.argmin(abs(ncdf.variables["nav_lon"][0,:]-(np.max(corners[0][:])+halo)))]
    jj=[np.argmin(abs(ncdf.variables["nav_lat"][:,0]-(np.min(corners[1][:])-halo))),
        np.argmin(abs(ncdf.variables["nav_lat"][:,0]-(np.max(corners[1][:])+halo)))]

    return min(ii), max(ii), min(jj), max(jj)
#
# Get the corners for a regular lat/lon grid
#
def corners(lon, lat) :
    jjm,iim = lon.shape
    cornerspoly = []
    cornersll = []
35 36
    centersll=[]
    radiusll = []
37
    index = []
38 39
    hdlon=np.mean(np.abs(np.diff(lon[0,:])))
    hdlat=np.mean(np.abs(np.diff(lat[:,0])))
40 41 42
    #
    for i in range(iim) :
        for j in range(jjm) :
43
            #
44
            boxll = RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat, 2)
45 46
            lons = [p[0] for p in boxll]
            lats = [p[1] for p in boxll]
47
            #
48
            cornerspoly.append(polygon.SphericalPolygon.from_lonlat(lons, lats, center=[lon[j,i], lat[j,i]]))
49
            #
50 51
            centersll.append([lon[j,i], lat[j,i]])
            radiusll.append(RPP.maxradius([lon[j,i], lat[j,i]], lons, lats))
52
            #
53
            index.append([j,i])
54
            cornersll.append(boxll)
55

56
    return cornersll, cornerspoly, centersll, radiusll, index
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
#
def gather(x, index) :
    y=[]
    for ia in index :
        y.append(list(x[ia[0][i],ia[1][i]] for i in range(len(ia[0]))))
    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 :
    def __init__(self, lolacorners) :
        #
        filename=config.get("OverAll", "HydroFile")
        INFO("Opening in HydroGrid :"+filename)
        self.ncfile=Dataset(filename,'r')
        istr, iend, jstr, jend = getbox(self.ncfile, lolacorners)
        self.box=[istr, iend, jstr, jend]
        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 87 88 89 90 91 92
        self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat)

    def select(self, c, r) :
        indices=[]
        for i in range(len(self.centers)) :
            if RPP.loladist(c,self.centers[i]) <= r+self.radius[i] :
                indices.append(i)
        return indices
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

class HydroData :
    def __init__(self, nf, box, index) :
        istr, iend, jstr, jend = box[:]
        self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index)
        self.tripdesc=nf.variables["trip"].long_name
        #
        self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index)
        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")
            self.basinsmax=np.max(np.where(nf.variables["basins"][:,:] <  1.e10))
        #
        self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index)
        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))
        #
        #
        self.topoindh=gather(nf.variables["topoind_h"][jstr:jend,istr:iend].astype(np.float32), index)
        self.topoindhdesc=nf.variables["topoind_h"].long_name
        att = getattrcontaining(nf, "topoind_h", "min")
        if len(att) > 0 :
            self.topoindhmin=att[0]
        else :
            INFO("We need to scan full file to find minimum topoind_h over domain")
            self.topoindhmin=np.min(np.where(nf.variables["topoind_h"][:,:] <  1.e15))
        #
        self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index)
        self.distodesc=nf.variables["disto"].long_name
        #
        self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index)
        self.facdesc=nf.variables["fac"].long_name
133 134 135 136 137 138 139 140 141 142 143 144
        #
        if "orog" in nf.variables.keys():
            self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index)
            self.orogdesc=nf.variables["orog"].long_name
        else:
            self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
        #
        if "floodplains" in nf.variables.keys():
            self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index)
            self.floodplainsdesc=nf.variables["floodplains"].long_name
        else:
            self.floodplains = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)