#
# 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 os
import sys
#
import getargs
config = getargs.SetupConfig()
EarthRadius=config.getfloat("OverAll", "EarthRadius", fallback=6370000.0)
FloodplainsFile=config.get("OverAll", "FloodplainsFile", fallback=None)
#
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) :
#
# Check that the domain fits in the data available.
lonrange=[getattrcontaining(ncdf, "nav_lon", "min"), getattrcontaining(ncdf, "nav_lon", "max")]
latrange=[getattrcontaining(ncdf, "nav_lat", "min"), getattrcontaining(ncdf, "nav_lat", "max")]
if len(lonrange[0]) <= 0 :
lonrange=[np.min(ncdf.variables["nav_lon"][:,:]),np.max(ncdf.variables["nav_lon"][:,:])]
if len(latrange[0]) <= 0 :
latrange=[np.min(ncdf.variables["nav_lat"][:,:]),np.max(ncdf.variables["nav_lat"][:,:])]
if np.min(corners[0][:]) < np.min(lonrange) or \
np.max(corners[0][:]) > np.max(lonrange) or \
np.min(corners[1][:]) < np.min(latrange) or \
np.max(corners[1][:]) > np.max(latrange) :
ERROR("The atmospheric domain does not fit inside the area covered by HydroFile. Please check the file.")
sys.exit()
#
# 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][:])))]
return min(ii)-halo_pts, max(ii)+halo_pts, min(jj)-halo_pts, max(jj)+halo_pts
#
# Get the corners for a regular lat/lon grid
#
def centers(lon, lat) :
jjm,iim = lon.shape
#
index = [[j,i] for i in range(iim) for j in range(jjm)]
centersll = [[lon[j,i], lat[j,i]] for j,i in index]
return index, centersll
def corners(subcentersll, hdlon, hdlat) :
cornersll = [RPP.boxit(np.array([lo, la]), hdlon, hdlat, 2) for lo,la in subcentersll]
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]
cornerspoly = [ polygon.SphericalPolygon.from_lonlat(lons, lats, center=cent) for lons, lats, cent in zip(Llons, Llats, subcentersll)]
#radiusll = [RPP.maxradius(np.array(cent), np.array(lons), np.array(lats)) for cent, lons, lats in zip(subcentersll,Llons,Llats)]
#
return cornersll, cornerspoly
#
def gather(x, index, default = 0) :
y=[]
for ia in index :
try:
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]) ))
except:
print("Error - the hydrological file may not recover all the domain", ia)
sys.exit()
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
#
def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) :
# The routing code in Fortran (i=1, j=2)
# 8 --- 1 --- 2
# | | |
# 7 --- X --- 3
# | | |
# 6 --- 5 --- 4
#
inc = np.array([[0,1,1,1,0,-1,-1,-1],[-1,-1,0,1,1,1,0,-1]])
elevation_precision=0.1
#
jlen,ilen = nf.variables["orog"].shape
if "missing_value" in nf.variables["orog"].ncattrs() :
missing = nf.variables["orog"].missing_value
elif "_FillValue" in nf.variables["orog"].ncattrs() :
missing = nf.variables["orog"]._FillValue
else :
print("No Missing values or Fill values in variable orography")
sys.exit()
land = np.where(nf.variables["orog"][jstr:jend,istr:iend] < missing)
topoind=np.zeros((jend-jstr,iend-istr))
topoind[:,:] = np.nan
#
for j,i in zip(land[:][0],land[:][1]) :
trip = int(nf.variables["trip"][jstr+j,istr+i]-1)
if trip < 8 :
inf=(istr+i+inc[0,trip])%ilen
jnf=(jstr+j+inc[1,trip])%jlen
dz=max(elevation_precision,nf.variables["orog"][jstr+j,istr+i]-nf.variables["orog"][jnf,inf])
dl=nf.variables["disto"][jstr+j,istr+i]-nf.variables["disto"][jnf,inf]
else :
if trip > 97 :
# Going into the ocean which is at elevation 0
dz = max(elevation_precision,nf.variables["orog"][jstr+j,istr+i])
else :
# Endorehic basin so we use thte lowest dz possible.
dz = elevation_precision
#
if nf.variables["disto"][jstr+j,istr+i] > 0 :
dl = nf.variables["disto"][jstr+j,istr+i]
else :
# Some other reasonable assumption
rx,ry = hydrogrid.resolution(i,j)
dl = (rx+ry)*0.5*0.5
topoind[j,i]=np.sqrt(dl**3./(dz*1.0e6))
#
minitopo = np.nanmin(topoind)
if minitopo <= np.finfo(float).eps :
ERROR("Topoind close to zero encoutered.")
sys.exit()
#
return gather(topoind, index, default = 10), "Topographic index which serves for the residence time", minitopo
#
class HydroGrid :
def __init__(self, lolacorners, wfile) :
#
self.source=config.get("OverAll", "HydroFile")
INFO("Opening in HydroGrid : "+self.source)
self.ncfile=Dataset(self.source,'r')
#
# Test if floodplains is correct
if FloodplainsFile is not None:
try:
nfp = Dataset(FloodplainsFile, "r")
njgfp, nigfp = nfp.variables["floodplains"].shape
njg, nig = self.ncfile.variables["nav_lon"].shape
if ((njg != njgfp) or (nig != nigfp)):
ERROR("Invalid FloodplainsFile format. Hydrogrod: {0} / FPfile: {1}".format((njg, nig),(njgfp, nigfp)))
sys.exit()
except:
ERROR("Invalid FloodplainsFile :{0}".format(FloodplainsFile))
sys.exit()
#
istr, iend, jstr, jend = getbox(self.ncfile, lolacorners)
self.hdlon = np.mean(np.abs(np.diff(self.ncfile.variables["nav_lon"][0,:])))
self.hdlat = np.mean(np.abs(np.diff(self.ncfile.variables["nav_lat"][:,0])))
self.hd = max(self.hdlon,self.hdlat)
# 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]
if jstr<0:
jstr = 0
if iend >=self.ncfile.variables["nav_lat"][0,:].shape[0]:
iend = self.ncfile.variables["nav_lat"][0,:].shape[0]
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("# Dimensions :"+str(self.iim)+" -- "+str(self.jjm))
if self.iim > 0 and self.jjm > 0 :
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)))
#
if wfile is None or not os.path.exists(wfile):
self.index, self.centers = centers(self.lon, self.lat)
def select(self, c, r) :
indices = [i for i in range(len(self.centers)) if (RPP.loladist(np.array(c),np.array(self.centers[i])) <= r+1.1*self.hd)]
subcenters = [self.centers[i] for i in indices]
polyll, polylist = corners(subcenters, self.hdlon, self.hdlat)
return indices, polylist
def resolution(self,i,j) :
if i+1 >= iim :
inx=i-1
else :
inx=i+1
if j+1 >= jjm :
jnx=j-1
else :
jnx=j+1
dlon = np.radians(np.abs(self.lon[j,i]-self.lon[j,inx]))
alon = np.cos(np.radians(self.lat[j,i]))**2*(np.sin(dlon/2))**2
clon = 2*np.arctan2(np.sqrt(alon), np.sqrt(1-alon))
dlat = np.radians(np.abs(self.lat[j,i]-self.lat[jnx,i]))
alat = (np.sin(dlat/2))**2
clat = 2*np.arctan2(np.sqrt(alat), np.sqrt(1-alat))
return EarthRadius*clon, EarthRadius*clat
class HydroData :
def __init__(self, nf, box, index, part, hydrogrid) :
istr, iend, jstr, jend = box[:]
#
# Flow direction
#
self.trip=gather(nf.variables["trip"][jstr:jend,istr:iend].astype(np.float32), index, default = 98)
self.tripdesc=nf.variables["trip"].long_name
#
# ID of basin
#
self.basins=gather(nf.variables["basins"][jstr:jend,istr:iend], index, 999)
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 the largest basin ID")
# This variable seems not to be used further
##self.basinsmax = part.domainmax(np.ma.max(np.ma.masked_where(self.basins < 1.e10, self.basins)))
##self.basinsmax = part.domainmax(np.max(np.array(self.basins)[np.array(self.basins) < missing]))
self.basinsmax = 999
#
# Distance to ocean
#
self.disto=gather(nf.variables["disto"][jstr:jend,istr:iend].astype(np.float32), index, default = 500)
self.distodesc=nf.variables["disto"].long_name
#
# Flow accumulation
#
self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index, default = 0)
self.facdesc=nf.variables["fac"].long_name
#
# 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.
#
if "orog" in nf.variables.keys() and "disto" in nf.variables.keys() :
self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index, default = 0)
self.orogdesc=nf.variables["orog"].long_name
self.topoind, self.topoinddesc, mintopo = calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid)
self.topoindmin = part.domainmin(mintopo)
elif "topoind" in nf.variables.keys() and "orog" in nf.variables.keys() :
self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, default = 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))
else :
ERROR("Not enough information in the HydroFile in order to compute topographic residence time.")
sys.exit()
#
#
if FloodplainsFile is not None:
nfp = Dataset(FloodplainsFile, "r")
self.floodplains = gather(nfp.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index, 0)
#self.floodplainsdesc=nfp.variables["floodplains"].long_name
elif "floodplains" in nf.variables.keys():
self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index, 0)
self.floodplainsdesc=nf.variables["floodplains"].long_name
else:
self.floodplains = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
class HydroParameter :
# Class to set, compute and manage the time constants of the routing scheme.
def __init__(self, hydrogrid) :
if hydrogrid.hd >= 0.5 :
# Case of the Fekete/Vörösmarty hydrological data set
self.stream_tcst = 0.24
self.fast_tcst = 3.0
self.slow_tcst = 25.0
self.flood_tcst = 4.0
self.swamp_cst = 0.2
elif hydrogrid.hd >= 0.016 :
# Case for MERIT
self.stream_tcst = 0.07
self.fast_tcst = 1.0
self.slow_tcst = 10.0
self.flood_tcst = 4.0
self.swamp_cst = 0.2
elif hydrogrid.hd >= 0.008 :
# Case for HydroSEHDS
self.stream_tcst = 0.01
self.fast_tcst = 0.5
self.slow_tcst = 7.0
self.flood_tcst = 4.0
self.swamp_cst = 0.2
else :
print("For the resolution ",hydrogrid.hd," deg. we do not yet have a parameter set.")
sys.exit