# -------------------------------------------------------------------------
# Name: READ METEO input maps
# Purpose:
#
# Author: PB
#
# Created: 13/07/2016
# Copyright: (c) PB 2016
# -------------------------------------------------------------------------
from cwatm.management_modules.data_handling import *
import scipy.ndimage
from scipy.interpolate import RegularGridInterpolator
[docs]class readmeteo(object):
"""
READ METEOROLOGICAL DATA
reads all meteorological data from netcdf4 files
**Global variables**
===================================== ====================================================================== =====
Variable [self.var] Description Unit
===================================== ====================================================================== =====
DtDay seconds in a timestep (default=86400) s
con_precipitation conversion factor for precipitation --
con_e conversion factor for evaporation --
ETRef potential evapotranspiration rate from reference crop m
Precipitation Precipitation (input for the model) m
only_radition --
TMin minimum air temperature K
TMax maximum air temperature K
Tavg Input, average air Temperature K
Rsds short wave downward surface radiation fluxes W/m2
EAct --
Psurf Instantaneous surface pressure Pa
Qair specific humidity kg/kg
Rsdl long wave downward surface radiation fluxes W/m2
Wind wind speed m/s
EWRef potential evaporation rate from water surface m
meteomapsscale if meteo maps have the same extend as the other spatial static maps -> --
meteodown if meteo maps should be downscaled --
InterpolationMethod --
buffer --
preMaps choose between steady state precipitation maps for steady state modflo --
tempMaps choose between steady state temperature maps for steady state modflow --
evaTMaps choose between steady state ETP water maps for steady state modflow or --
eva0Maps choose between steady state ETP reference maps for steady state modflo --
glaciermeltMaps --
glacierrainMaps --
wc2_tavg High resolution WorldClim map for average temperature K
wc4_tavg upscaled to low resolution WorldClim map for average temperature K
wc2_tmin High resolution WorldClim map for min temperature K
wc4_tmin upscaled to low resolution WorldClim map for min temperature K
wc2_tmax High resolution WorldClim map for max temperature K
wc4_tmax upscaled to low resolution WorldClim map for max temperature K
wc2_prec High resolution WorldClim map for precipitation m
wc4_prec upscaled to low resolution WorldClim map for precipitation m
xcoarse_prec --
ycoarse_prec --
xfine_prec --
yfine_prec --
meshlist_prec --
xcoarse_tavg --
ycoarse_tavg --
xfine_tavg --
yfine_tavg --
meshlist_tavg --
meteo --
prec precipitation in m m
temp average temperature in Celsius deg °C
WtoMJ Conversion factor from [W] to [MJ] for radiation: 86400 * 1E-6 --
includeGlaciers --
includeOnlyGlaciersMelt --
GlacierMelt --
GlacierRain --
===================================== ====================================================================== =====
**Functions**
"""
def __init__(self, model):
self.model = model
self.var = model.var
[docs] def initial(self):
"""
Initial part of meteo
read multiple file of input
"""
# fit meteorological forcing data to size and resolution of mask map
#-------------------------------------------------------------------
name = cbinding('PrecipitationMaps')
nameall = glob.glob(os.path.normpath(name))
if not nameall:
msg = "Error 215: In readmeteo, cannot find precipitation maps "
raise CWATMFileError(name, msg, sname='PrecipitationMaps')
namemeteo = nameall[0]
latmeteo, lonmeteo, cell, invcellmeteo, rows, cols = readCoordNetCDF(namemeteo)
nameldd = cbinding('Ldd')
#nameldd = os.path.splitext(nameldd)[0] + '.nc'
#latldd, lonldd, cell, invcellldd, row, cols = readCoordNetCDF(nameldd)
latldd, lonldd, cell, invcellldd, rows, cols = readCoord(nameldd)
maskmapAttr['reso_mask_meteo'] = round(invcellldd / invcellmeteo)
# if meteo maps have the same extend as the other spatial static maps -> meteomapsscale = True
self.var.meteomapsscale = True
if invcellmeteo != invcellldd:
if (not(Flags['quiet'])) and (not(Flags['veryquiet'])) and (not(Flags['check'])):
msg = "Resolution of meteo forcing is " + str(maskmapAttr['reso_mask_meteo']) + " times higher than base maps."
print(msg)
self.var.meteomapsscale = False
cutmap[0], cutmap[1], cutmap[2], cutmap[3] = mapattrNetCDF(nameldd)
for i in range(4): cutmapFine[i] = cutmap[i]
# for downscaling meteomaps , Wordclim data at a finer resolution is used
# here it is necessary to clip the wordclim data so that they fit to meteo dataset
self.var.meteodown = False
# if interpolationmethod not defined in settingsfil, use spline interpolation
self.var.InterpolationMethod = 'spline'
self.var.buffer = False
if "usemeteodownscaling" in binding:
self.var.meteodown = returnBool('usemeteodownscaling')
if 'InterpolationMethod' in binding:
# interpolation option can be spline or bilinear
self.var.InterpolationMethod = cbinding('InterpolationMethod')
if self.var.InterpolationMethod != 'bilinear' and self.var.InterpolationMethod != 'spline' and self.var.InterpolationMethod != 'kron':
msg = 'Error: InterpolationMethod in settings file must be one of the following: "spline" or "bilinear", but it is {}'.format(self.var.InterpolationMethod)
raise CWATMError(msg)
if self.var.InterpolationMethod == 'bilinear':
self.var.buffer = True
check_clim = False
if self.var.meteodown:
check_clim = checkMeteo_Wordclim(namemeteo, cbinding('downscale_wordclim_prec'))
# in case other mapsets are used e.g. Cordex RCM meteo data
if (latldd != latmeteo) or (lonldd != lonmeteo):
cutmapFine[0], cutmapFine[1], cutmapFine[2], cutmapFine[3], cutmapVfine[0], cutmapVfine[1], cutmapVfine[2], cutmapVfine[3] = mapattrNetCDFMeteo(namemeteo)
if not self.var.meteomapsscale:
# if the cellsize of the spatial dataset e.g. ldd, soil etc is not the same as the meteo maps than:
cutmapFine[0], cutmapFine[1],cutmapFine[2],cutmapFine[3],cutmapVfine[0], cutmapVfine[1],cutmapVfine[2],cutmapVfine[3] = mapattrNetCDFMeteo(namemeteo)
# downscaling wordlclim maps
for i in range(4): cutmapGlobal[i] = cutmapFine[i]
if not(check_clim):
# for downscaling it is always cut from the global map
if (latldd != latmeteo) or (lonldd != lonmeteo):
cutmapGlobal[0] = int(cutmap[0] / maskmapAttr['reso_mask_meteo'])
cutmapGlobal[2] = int(cutmap[2] / maskmapAttr['reso_mask_meteo'])
cutmapGlobal[1] = int(cutmap[1] / maskmapAttr['reso_mask_meteo']+0.999)
cutmapGlobal[3] = int(cutmap[3] / maskmapAttr['reso_mask_meteo']+0.999)
# -------------------------------------------------------------------
self.var.includeGlaciers = False
self.var.includeOnlyGlaciersMelt = False
if 'includeGlaciers' in option:
self.var.includeGlaciers = checkOption('includeGlaciers')
if 'includeOnlyGlaciersMelt' in option:
self.var.includeOnlyGlaciersMelt = checkOption('includeOnlyGlaciersMelt')
self.var.preMaps = 'PrecipitationMaps'
self.var.tempMaps = 'TavgMaps'
self.var.evaTMaps = 'ETMaps'
self.var.eva0Maps = 'E0Maps'
if self.var.includeGlaciers:
self.var.glaciermeltMaps = 'MeltGlacierMaps'
if not self.var.includeOnlyGlaciersMelt:
self.var.glacierrainMaps = 'PrecGlacierMaps'
self.var.only_radition = False
if 'only_radiation' in binding:
self.var.only_radition = returnBool('only_radiation')
if checkOption('calc_evaporation'):
if self.var.only_radition:
# for maps from EMO-5 with total radiation and vapor pressure instead of huss, air pressure, rsds and rlds
meteomaps = [self.var.preMaps, self.var.tempMaps,'TminMaps','TmaxMaps','WindMaps','RGDMaps','EActMaps']
else:
meteomaps = [self.var.preMaps, self.var.tempMaps,'TminMaps','TmaxMaps','PSurfMaps','WindMaps','RSDSMaps','RSDLMaps']
if returnBool('useHuss'):
meteomaps.append('QAirMaps')
else:
meteomaps.append('RhsMaps')
if self.var.includeGlaciers:
meteomaps.append(self.var.glaciermeltMaps)
if not self.var.includeOnlyGlaciersMelt:
meteomaps.append(self.var.glacierrainMaps)
else:
if self.var.includeGlaciers:
if not self.var.includeOnlyGlaciersMelt:
meteomaps = [self.var.preMaps, self.var.tempMaps, self.var.evaTMaps, self.var.eva0Maps, self.var.glaciermeltMaps, self.var.glacierrainMaps]
else:
meteomaps = [self.var.preMaps, self.var.tempMaps, self.var.evaTMaps, self.var.eva0Maps,
self.var.glaciermeltMaps]
else:
meteomaps = [self.var.preMaps, self.var.tempMaps,self.var.evaTMaps,self.var.eva0Maps]
#meteomaps = ["PrecipitationMaps","TavgMaps"]
multinetdf(meteomaps)
# downscaling to wordclim, set parameter to 0 in case they are only used as dummy
self.var.wc2_tavg = 0
self.var.wc4_tavg = 0
self.var.wc2_tmin = 0
self.var.wc4_tmin = 0
self.var.wc2_tmax = 0
self.var.wc4_tmax = 0
self.var.wc2_prec = 0
self.var.wc4_prec = 0
if self.var.InterpolationMethod == 'bilinear':
#these variables are generated to avoid calculating them at each timestep
self.var.xcoarse_prec = 0
self.var.ycoarse_prec = 0
self.var.xfine_prec = 0
self.var.yfine_prec = 0
self.var.meshlist_prec = 0
self.var.xcoarse_tavg = 0
self.var.ycoarse_tavg = 0
self.var.xfine_tavg = 0
self.var.yfine_tavg = 0
self.var.meshlist_tavg = 0
# read dem for making a anomolydem between high resolution dem and low resoultion dem
"""
# for downscaling1
dem = loadmap('Elevation', compress = False, cut = False)
demHigh = dem[cutmapFine[2]*6:cutmapFine[3]*6, cutmapFine[0]*6:cutmapFine[1]*6]
rows = demHigh.shape[0]
cols = demHigh.shape[1]
dem2 = demHigh.reshape(rows/6,6,cols/6,6)
dem3 = np.average(dem2, axis=(1, 3))
demLow = np.kron(dem3, np.ones((6, 6)))
demAnomaly = demHigh - demLow
self.var.demHigh = compressArray(demHigh[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]],pcr = False)
self.var.demAnomaly = compressArray(demAnomaly[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]],pcr = False)
"""
# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
#def downscaling1(self,input, downscale = 0):
"""
Downscaling based on elevation correction for temperature and pressure
:param input:
:param downscale: 0 for no change, 1: for temperature change 6 deg per 1km , 2 for psurf
:return: input - downscaled input data
"""
"""
# if meteo maps have the same extend as the other spatial static maps -> meteomapsscale = True
if not self.var.meteomapsscale:
down1 = np.kron(input, np.ones((6, 6)))
down2 = down1[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]].astype(np.float64)
down3 = compressArray(down2)
if downscale == 0:
input = down3
if downscale == 1:
# temperature scaling 6 deg per 1000m difference in altitude
# see overview in Minder et al 2010 - http://onlinelibrary.wiley.com/doi/10.1029/2009JD013493/full
tempdiff = -0.006 * self.var.demAnomaly
input = down3 + tempdiff
if downscale == 2:
# psurf correction
# https://www.sandhurstweather.org.uk/barometric.pdf
# factor = exp(-elevation / (Temp x 29.263) Temp in deg K
demLow = self.var.demHigh - self.var.demAnomaly
tavgK = self.var.Tavg + 273.15
factor1 = np.exp(-1 * demLow / (tavgK * 29.263))
factor2 = np.exp(-1 * self.var.demHigh / (tavgK * 29.263))
sealevelpressure = down3 / factor1
input = sealevelpressure * factor2
return input
"""
# def downscaling2_peter(self,input, downscaleName = "", wc2 = 0 , wc4 = 0, x=None, y=None, xfine=None, yfine=None, meshlist=None, downscale = 0):
# """
# Downscaling with only internal (inside the coarse gridcell) interpolation
#
# :param input: low input map
# :param downscaleName: High resolution monthly map from WorldClim
# :param wc2: High resolution WorldClim map
# :param wc4: upscaled to low resolution
# :param downscale: 0 for no change, 1: for temperature , 2 for pprecipitation, 3 for psurf
# :return: input - downscaled input data
# :return: wc2
# :return: wc4
# """
# reso = maskmapAttr['reso_mask_meteo']
# resoint = int(reso)
# if self.var.meteomapsscale:
# if downscale == 0:
# return input
# else:
# return input, wc2, wc4
#
# down3 = np.kron(input, np.ones((resoint, resoint)))
# # this is creating an array resoint times bigger than input, by copying each item resoint times in x and y direction
#
# if downscale == 0:
# down2 = down3[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]].astype(np.float64)
# input = compressArray(down2)
# return input
# else:
# if dateVar['newStart'] or dateVar['newMonth']: # loading every month a new map
# # wc1 = readnetcdf2(downscaleName, dateVar['currDate'], useDaily='month', compress = False, cut = False)
# # wc2 = wc1[cutmapGlobal[2]*resoint:cutmapGlobal[3]*resoint, cutmapGlobal[0]*resoint:cutmapGlobal[1]*resoint]
# #print('\n'.join([' '.join(['{:4}'.format(item) for item in row]) for row in wc2]))
#
# if downscale == 2: # precipitation
# #wc3 looks a like wc3
# wc3 = wc2.reshape(wc2.shape[0] // resoint, resoint, wc2.shape[1] // resoint, resoint)
# #wc3mean looks like w4
# wc3mean = np.nanmean(wc3, axis=(1, 3))
# # Average of wordclim on the bigger input raster scale
# wc3kron = np.kron(wc3mean, np.ones((resoint, resoint)))
# # the average values are spread out to the fine scale
# #looks like quot_wc, but wc2 = input, wc3kron = wc4
# wc4 = divideValues(wc2, wc3kron)
# # wc4 holds the correction multiplicator on fine scale
#
# if downscale == 1: # Temperature
# #diff wc is different because originally it is wc4 - input
# #diff_wc is difference on small scale
# diff_wc = wc2 - down3
# # on fine scale: wordclim fine scale - spreaded input data (same value for each big cell)
# wc3 = diff_wc.reshape(wc2.shape[0] // resoint, resoint, wc2.shape[1] // resoint, resoint)
# wc4 = np.nanmean(wc3, axis=(1, 3))
# wc4kron = np.kron(wc4, np.ones((resoint, resoint)))
# # wordclim is averaged on big cell scale and the average is spread out to fine raster
# down1 = diff_wc - wc4kron + down3
# # result is the fine scale input data + the difference of wordclim - input data - the average difference of wordclim - input
# down1 = np.where(np.isnan(down1),down3,down1)
# if downscale == 2: # precipitation
# # in the other interpolations this is wc2 * quotSmooth, wc2 being the fine worldclimmap cut to map extent, quotSmooth being the interpolated difference between the input and summed worldclim
# down1 = down3 * wc4
# down1 = np.where(np.isnan(down1),down3,down1)
# down1 = np.where(np.isinf(down1), down3, down1)
#
# down2 = down1[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]].astype(np.float64)
# input = compressArray(down2)
# return input, wc2, wc4
# --- end downscaling ----------------------------
[docs] def downscaling2(self,input, downscaleName = "", wc2 = 0 , wc4 = 0, x=None, y=None, xfine=None, yfine=None, meshlist=None, MaskMapBoundaries= None, downscale = 0):
"""
Downscaling based on Delta method:
Note:
| **References**
| Moreno and Hasenauer 2015:
| ftp://palantir.boku.ac.at/Public/ClimateData/Moreno_et_al-2015-International_Journal_of_Climatology.pdf
| Mosier et al. 2018:
| http://onlinelibrary.wiley.com/doi/10.1002/joc.5213/epdf\
:param input: low input map
:param downscaleName: High resolution monthly map from WorldClim
:param wc2: High resolution WorldClim map
:param wc4: upscaled to low resolution
:param MaskMapBoundaries: if 1 maskmap does not touch meteo input dataset boundary, if 0 maskmap touches it
:param downscale: 0 for no change, 1: for temperature , 2 for pprecipitation, 3 for psurf
:return: input - downscaled input data
:return: wc2
:return: wc4
"""
reso = maskmapAttr['reso_mask_meteo']
resoint = int(reso)
if self.var.InterpolationMethod == 'bilinear' and (downscale == 1 or downscale == 2):
buffer = 1
buffer1, buffer2, buffer3, buffer4 = MaskMapBoundaries
#if 1: does not touch boundaries of meteo input map, if 0 touches boundary of input map
# buffer1
# ---------
# buffer3¦ ¦ buffer4
# ¦ ¦
# ---------
# buffer2
# to perform bilinear interpolation a buffer around the maskmap is needed, if maskmap touches bounary of input map an artifical buffer has to be created by duplicating the last row/column
if buffer1 == 0:
input_first_row = input[0, :]
input = np.vstack((input_first_row[np.newaxis, :], input))
if buffer2 == 0:
input_last_row = input[-1, :]
input = np.vstack((input, input_last_row[np.newaxis, :]))
if buffer3 == 0:
input_first_column = input[:, 0]
input = np.hstack((input_first_column[:, np.newaxis], input))
if buffer4 == 0:
input_last_column = input[:, -1]
input = np.hstack((input, input_last_column[:, np.newaxis]))
if dateVar['newStart']:
x = np.arange(0.5, np.shape(input)[0] + 0.5)
y = np.arange(0.5, np.shape(input)[1] + 0.5)
xfine = np.arange(0.5 + 1 / (resoint * 2), np.shape(input)[0] - 0.5, 1 / resoint)
yfine = np.arange(0.5 + 1 / (resoint * 2), np.shape(input)[1] - 0.5, 1 / resoint)
xmesh, ymesh = np.meshgrid(xfine, yfine)
meshlist = list(zip(xmesh.flatten(), ymesh.flatten()))
else:
buffer = 0
if self.var.meteomapsscale:
if downscale == 0:
return input
else:
return input, wc2, wc4
if buffer == 0:
# this is creating an array resoint times bigger than input, by copying each item resoint times in x and y direction
down3 = np.kron(input, np.ones((resoint, resoint)))
else:
down3 = np.kron(input[buffer:-buffer, buffer:-buffer], np.ones((resoint, resoint)))
if downscale == 0:
down2 = down3[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]].astype(np.float64)
input = compressArray(down2)
return input
else:
if dateVar['newStart'] or dateVar['newMonth']: # loading every month a new map
wc1 = readnetcdf2(downscaleName, dateVar['currDate'], useDaily='month', compress = False, cut = False)
if self.var.InterpolationMethod == 'bilinear':
# to perform bilinear interpolation a buffer around the maskmap is needed, if maskmap touches bounary of input map an artifical buffer has to be created by duplicating the last row/column
if buffer1 == 0:
wc1_first_row = wc1[:resoint, :]
wc1 = np.vstack((wc1_first_row, wc1))
if buffer2 == 0:
wc1_last_row = wc1[-resoint:, :]
wc1 = np.vstack((wc1, wc1_last_row))
if buffer3 == 0:
wc1_first_column = wc1[:, :resoint]
wc1 = np.hstack((wc1_first_column, wc1))
if buffer4 == 0:
wc1_last_column = wc1[:, -resoint:]
wc1 = np.hstack((wc1, wc1_last_column))
#include buffer
if buffer1 == 0:
#if maskmap reaches upper boundary you cannot do -buffer because than index would be negative
wc2 = wc1[int(np.floor((cutmapGlobal[2]) * reso)):int(np.ceil((cutmapGlobal[3] + buffer* 2) * reso)), :]
if buffer3 == 0:
wc2 = wc2[:, int(np.floor((cutmapGlobal[0]) * reso)): int(
np.ceil((cutmapGlobal[1] + buffer * 2) * reso))]
elif buffer3 == 0:
# if maskmap reaches left boundary you cannot do -buffer because than index would be negative
wc2 = wc1[int(np.floor((cutmapGlobal[2] - buffer) * reso)):int(
np.ceil((cutmapGlobal[3] + buffer) * reso)),int(np.floor((cutmapGlobal[0]) * reso)) : int(np.ceil((cutmapGlobal[1] + buffer * 2) * reso))]
else:
wc2 = wc1[int(np.floor((cutmapGlobal[2] - buffer) * reso)):int(
np.ceil((cutmapGlobal[3] + buffer) * reso)),
int(np.floor((cutmapGlobal[0] - buffer) * reso)):int(
np.ceil((cutmapGlobal[1] + buffer) * reso))]
else:
wc2 = wc1[(cutmapGlobal[2] - buffer) * resoint: (cutmapGlobal[3] + buffer) * resoint,
(cutmapGlobal[0] - buffer) * resoint: (cutmapGlobal[1] + buffer) * resoint]
rows = wc2.shape[0]
cols = wc2.shape[1]
wc3 = wc2.reshape(rows//resoint,resoint,cols//resoint,resoint)
wc4 = np.nanmean(wc3, axis=(1, 3))
# wc4 is as big as the input array -> average of the fine scale downscale map
if self.var.InterpolationMethod == 'kron':
if downscale == 2: # precipitation
# wc3 looks a like wc3
wc3 = wc2.reshape(wc2.shape[0] // resoint, resoint, wc2.shape[1] // resoint, resoint)
# wc3mean looks like w4
wc3mean = np.nanmean(wc3, axis=(1, 3))
# Average of wordclim on the bigger input raster scale
wc3kron = np.kron(wc3mean, np.ones((resoint, resoint)))
# the average values are spread out to the fine scale
# looks like quot_wc, but wc2 = input, wc3kron = wc4
wc4 = divideValues(wc2, wc3kron)
if downscale == 1: # Temperature
diff_wc = wc4 - input
if self.var.InterpolationMethod == 'spline':
diffSmooth = scipy.ndimage.zoom(diff_wc, resoint, order=1)
down1 = wc2 - diffSmooth
elif self.var.InterpolationMethod == 'bilinear':
bilinear_interpolation = RegularGridInterpolator((x, y), diff_wc)
diffSmooth = bilinear_interpolation(meshlist)
diffSmooth = diffSmooth.reshape(len(xfine), len(yfine), order='F')
#no buffer for real downscaled values
crop = int(resoint / 2)
diffSmooth = diffSmooth[crop:-crop, crop:-crop]
down1 = wc2[buffer * resoint:-buffer * resoint, buffer * resoint:-buffer * resoint] - diffSmooth
elif self.var.InterpolationMethod == 'kron':
diff_wc = wc2 - down3
# on fine scale: wordclim fine scale - spreaded input data (same value for each big cell)
wc3 = diff_wc.reshape(wc2.shape[0] // resoint, resoint, wc2.shape[1] // resoint, resoint)
wc4 = np.nanmean(wc3, axis=(1, 3))
wc4kron = np.kron(wc4, np.ones((resoint, resoint)))
# wordclim is averaged on big cell scale and the average is spread out to fine raster
down1 = diff_wc - wc4kron + down3
# result is the fine scale input data + the difference of wordclim - input data - the average difference of wordclim - input
#down1 = np.where(np.isnan(down1), down3, down1)
down1 = np.where(np.isnan(down1),down3,down1)
if downscale == 2: # precipitation
if self.var.InterpolationMethod == 'spline':
quot_wc = divideValues(input, wc4)
quotSmooth = scipy.ndimage.zoom(quot_wc, resoint, order=1)
down1 = wc2 * quotSmooth
elif self.var.InterpolationMethod == 'bilinear':
quot_wc = divideValues(input, wc4)
bilinear_interpolation = RegularGridInterpolator((x, y), quot_wc)
quotSmooth = bilinear_interpolation(meshlist)
quotSmooth = quotSmooth.reshape(len(xfine), len(yfine), order='F')
crop = int(resoint/2)
quotSmooth = quotSmooth[crop:-crop, crop:-crop]
down1 = wc2[buffer * resoint:-buffer * resoint, buffer * resoint:-buffer * resoint] * quotSmooth
elif self.var.InterpolationMethod == 'kron':
down1 = down3 * wc4
down1 = np.where(np.isnan(down1),down3,down1)
down1 = np.where(np.isinf(down1), down3, down1)
down2 = down1[cutmapVfine[2]:cutmapVfine[3], cutmapVfine[0]:cutmapVfine[1]].astype(np.float64)
input = compressArray(down2)
if self.var.InterpolationMethod == 'bilinear' and (downscale == 1 or downscale == 2):
return input, wc2, wc4, x, y, xfine, yfine, meshlist
return input, wc2, wc4
# --- end downscaling ----------------------------
[docs] def dynamic(self):
"""
Dynamic part of the readmeteo module
Read meteo input maps from netcdf files
Note:
If option *calc_evaporation* is False only precipitation, avg. temp., and 2 evaporation vlaues are read
Otherwise all the variable needed for Penman-Monteith
Note:
If option *TemperatureInKelvin* = True temperature is assumed to be Kelvin instead of Celsius!
"""
if Flags['warm']:
# if warmstart use stored meteo variables
no = dateVar['curr']-1
self.var.Precipitation = self.var.meteo[0,no]
self.var.Tavg = self.var.meteo[1,no]
self.var.ETRef = self.var.meteo[2,no]
self.var.EWRef = self.var.meteo[3,no]
if self.var.includeGlaciers:
self.var.GlacierMelt = self.var.meteo[4, no]
if not self.var.includeOnlyGlaciersMelt:
self.var.GlacierRain = self.var.meteo[5, no]
return
ZeroKelvin = 0.0
if checkOption('TemperatureInKelvin'):
# if temperature is in Kelvin -> conversion to deg C
# TODO in initial there could be a check if temperature > 200 -> automatic change to Kelvin
ZeroKelvin = 273.15
self.var.Precipitation, MaskMapBoundary = readmeteodata(self.var.preMaps, dateVar['currDate'], addZeros=True, mapsscale = self.var.meteomapsscale, buffering= self.var.buffer)
self.var.Precipitation = self.var.Precipitation * self.var.DtDay * self.var.con_precipitation
self.var.Precipitation = np.maximum(0., self.var.Precipitation)
if self.var.includeGlaciers:
self.var.GlacierMelt, MaskMapBoundary = readmeteodata(self.var.glaciermeltMaps, dateVar['currDate'], addZeros=True, mapsscale = True)
if not self.var.includeOnlyGlaciersMelt:
self.var.GlacierRain, MaskMapBoundary = readmeteodata(self.var.glacierrainMaps, dateVar['currDate'], addZeros=True, mapsscale = True)
if self.var.meteodown:
if self.var.InterpolationMethod == 'bilinear':
self.var.Precipitation, self.var.wc2_prec, self.var.wc4_prec, self.var.xcoarse_prec, self.var.ycoarse_prec, self.var.xfine_prec, self.var.yfine_prec, self.var.meshlist_prec = self.downscaling2(
self.var.Precipitation, "downscale_wordclim_prec", self.var.wc2_prec, self.var.wc4_prec,
self.var.xcoarse_prec, self.var.ycoarse_prec, self.var.xfine_prec, self.var.yfine_prec,
self.var.meshlist_prec, MaskMapBoundary, downscale=2)
else:
self.var.Precipitation, self.var.wc2_prec, self.var.wc4_prec = self.downscaling2(self.var.Precipitation, "downscale_wordclim_prec", self.var.wc2_prec, self.var.wc4_prec, downscale=2)
else:
self.var.Precipitation = self.downscaling2(self.var.Precipitation, "downscale_wordclim_prec", self.var.wc2_prec, self.var.wc4_prec, downscale=0)
#self.var.Precipitation = self.var.Precipitation * 1000
self.var.prec = self.var.Precipitation / self.var.con_precipitation
# precipitation (conversion to [m] per time step) `
if Flags['check']:
checkmap(self.var.preMaps, "", self.var.Precipitation, True, True, self.var.Precipitation)
#self.var.Tavg = readnetcdf2('TavgMaps', dateVar['currDate'], addZeros = True, zeros = ZeroKelvin, meteo = True)
tzero = 0
if checkOption('TemperatureInKelvin'):
tzero = ZeroKelvin
self.var.Tavg, MaskMapBoundary = readmeteodata(self.var.tempMaps,dateVar['currDate'], addZeros=True, zeros = tzero, mapsscale = self.var.meteomapsscale, buffering= self.var.buffer)
if self.var.meteodown:
if self.var.InterpolationMethod == 'bilinear':
self.var.Tavg, self.var.wc2_tavg, self.var.wc4_tavg, self.var.xcoarse_tavg, self.var.ycoarse_tavg, self.var.xfine_tavg, self.var.yfine_tavg, self.var.meshlist_tavg = self.downscaling2(
self.var.Tavg, "downscale_wordclim_tavg", self.var.wc2_tavg, self.var.wc4_tavg, self.var.xcoarse_tavg,
self.var.ycoarse_tavg, self.var.xfine_tavg, self.var.yfine_tavg, self.var.meshlist_tavg, MaskMapBoundary, downscale=1)
else:
self.var.Tavg, self.var.wc2_tavg, self.var.wc4_tavg = self.downscaling2(self.var.Tavg, "downscale_wordclim_tavg", self.var.wc2_tavg, self.var.wc4_tavg, downscale=1)
else:
self.var.Tavg = self.downscaling2(self.var.Tavg, "downscale_wordclim_tavg", self.var.wc2_tavg, self.var.wc4_tavg, downscale=0)
self.var.temp = self.var.Tavg.copy()
# average DAILY temperature (even if you are running the model
# on say an hourly time step) [degrees C]
if checkOption('TemperatureInKelvin'):
self.var.Tavg -= ZeroKelvin
if Flags['check']:
checkmap(self.var.tempMaps, "", self.var.Tavg, True, True, self.var.Tavg)
#self.var.Tavg = downscaling(self.var.Tavg, downscale = 0)
# -----------------------------------------------------------------------
# if evaporation has to be calculated load all the meteo map sets
# Temparture min, max; Windspeed, specific humidity or relative humidity
# psurf, radiation
# -----------------------------------------------------------------------
if checkOption('calc_evaporation'):
#self.var.TMin = readnetcdf2('TminMaps', dateVar['currDate'], addZeros = True, zeros = ZeroKelvin, meteo = True)
self.var.TMin, MaskMapBoundary = readmeteodata('TminMaps',dateVar['currDate'], addZeros=True, zeros=ZeroKelvin, mapsscale = self.var.meteomapsscale, buffering= self.var.buffer)
if self.var.meteodown:
if self.var.InterpolationMethod == 'bilinear':
self.var.TMin, self.var.wc2_tmin, self.var.wc4_tmin, _, _, _, _, _ = self.downscaling2(self.var.TMin,
"downscale_wordclim_tmin",
self.var.wc2_tmin,
self.var.wc4_tmin,
self.var.xcoarse_tavg,
self.var.ycoarse_tavg,
self.var.xfine_tavg,
self.var.yfine_tavg,
self.var.meshlist_tavg, MaskMapBoundary,
downscale=1)
else:
self.var.TMin, self.var.wc2_tmin, self.var.wc4_tmin = self.downscaling2(self.var.TMin, "downscale_wordclim_tmin", self.var.wc2_tmin, self.var.wc4_tmin, downscale=1)
else:
self.var.TMin = self.downscaling2(self.var.TMin, "downscale_wordclim_tmin", self.var.wc2_tmin, self.var.wc4_tmin, downscale=0)
if Flags['check']:
checkmap('TminMaps', "", self.var.TMin, True, True, self.var.TMin)
#self.var.TMax = readnetcdf2('TmaxMaps', dateVar['currDate'], addZeros = True, zeros = ZeroKelvin, meteo = True)
self.var.TMax, MaskMapBoundary = readmeteodata('TmaxMaps', dateVar['currDate'], addZeros=True, zeros=ZeroKelvin, mapsscale = self.var.meteomapsscale, buffering= self.var.buffer)
if self.var.meteodown:
if self.var.InterpolationMethod == 'bilinear':
self.var.TMax, self.var.wc2_tmax, self.var.wc4_tmax, _, _, _, _, _ = self.downscaling2(self.var.TMax,
"downscale_wordclim_tmin",
self.var.wc2_tmax,
self.var.wc4_tmax,
self.var.xcoarse_tavg,
self.var.ycoarse_tavg,
self.var.xfine_tavg,
self.var.yfine_tavg,
self.var.meshlist_tavg, MaskMapBoundary,
downscale=1)
else:
self.var.TMax, self.var.wc2_tmax, self.var.wc4_tmax = self.downscaling2(self.var.TMax, "downscale_wordclim_tmin", self.var.wc2_tmax, self.var.wc4_tmax, downscale=1)
else:
self.var.TMax = self.downscaling2(self.var.TMax, "downscale_wordclim_tmin", self.var.wc2_tmax, self.var.wc4_tmax, downscale=0)
if Flags['check']: checkmap('TmaxMaps', "", self.var.TMax, True, True, self.var.TMax)
if checkOption('TemperatureInKelvin'):
self.var.TMin -= ZeroKelvin
self.var.TMax -= ZeroKelvin
#self.var.Wind = readnetcdf2('WindMaps', dateVar['currDate'], addZeros = True, meteo = True)
self.var.Wind, MaskMapBoundary = readmeteodata('WindMaps', dateVar['currDate'], addZeros=True, mapsscale = self.var.meteomapsscale)
self.var.Wind = self.downscaling2(self.var.Wind)
# wind speed maps at 10m [m/s]
# Adjust wind speed for measurement height: wind speed measured at
# 10 m, but needed at 2 m height
# Shuttleworth, W.J. (1993) in Maidment, D.R. (1993), p. 4.36
self.var.Wind = self.var.Wind * 0.749
if self.var.only_radition:
# read daily calculated radiation [in KJ/m2/day]
# named here Rsds instead of rds, because use in evaproationPot in the same way as rsds
self.var.Rsds, MaskMapBoundary = readmeteodata('RGDMaps', dateVar['currDate'], addZeros=True, mapsscale=self.var.meteomapsscale)
#self.var.Rsds = self.downscaling2(self.var.Rsds) * 0.001 # convert from KJ to MJ/m2/day
self.var.Rsds = self.downscaling2(self.var.Rsds) * 0.000001 # convert from KJ to MJ/m2/day
# but for EMO it is 1e6 instead 1000 it seems it is J instead of KJ
# read daily vapor pressure [in hPa]
self.var.EAct, MaskMapBoundary = readmeteodata('EActMaps', dateVar['currDate'], addZeros=True, mapsscale=self.var.meteomapsscale)
self.var.EAct = self.downscaling2(self.var.EAct) * 0.1 # convert from hP to kP
return
#self.var.Psurf = readnetcdf2('PSurfMaps', dateVar['currDate'], addZeros = True, meteo = True)
self.var.Psurf, MaskMapBoundary = readmeteodata('PSurfMaps', dateVar['currDate'], addZeros=True, mapsscale = self.var.meteomapsscale)
self.var.Psurf = self.downscaling2(self.var.Psurf)
# Instantaneous surface pressure[Pa]
if returnBool('useHuss'):
#self.var.Qair = readnetcdf2('QAirMaps', dateVar['currDate'], addZeros = True, meteo = True)
self.var.Qair, MaskMapBoundary = readmeteodata('QAirMaps', dateVar['currDate'], addZeros=True, mapsscale =self.var.meteomapsscale)
# 2 m istantaneous specific humidity[kg / kg]
else:
#self.var.Qair = readnetcdf2('RhsMaps', dateVar['currDate'], addZeros = True, meteo = True)
self.var.Qair, MaskMapBoundary = readmeteodata('RhsMaps', dateVar['currDate'], addZeros=True, mapsscale =self.var.meteomapsscale)
self.var.Qair = self.downscaling2(self.var.Qair)
#self.var.Rsds = readnetcdf2('RSDSMaps', dateVar['currDate'], addZeros = True, meteo = True)
self.var.Rsds, MaskMapBoundary = readmeteodata('RSDSMaps', dateVar['currDate'], addZeros=True, mapsscale = self.var.meteomapsscale)
self.var.Rsds = self.downscaling2(self.var.Rsds)
# radiation surface downwelling shortwave maps [W/m2]
#self.var.Rsdl = readnetcdf2('RSDLMaps', dateVar['currDate'], addZeros = True, meteo = True)
self.var.Rsdl, MaskMapBoundary = readmeteodata('RSDLMaps', dateVar['currDate'], addZeros=True, mapsscale = self.var.meteomapsscale)
self.var.Rsdl = self.downscaling2(self.var.Rsdl)
# radiation surface downwelling longwave maps [W/m2]
#--------------------------------------------------------
# conversions
# [Pa] to [KPa]
self.var.Psurf = self.var.Psurf * 0.001
# Conversion factor from [W] to [MJ]
self.var.WtoMJ = 86400 * 1E-6
# conversion from W/m2 to MJ/m2/day
self.var.Rsds = self.var.Rsds * self.var.WtoMJ
self.var.Rsdl = self.var.Rsdl * self.var.WtoMJ
# if pot evaporation is already precalulated
else:
# in case ET_ref is the same resolution as the other meteo input map, there is an optional flag in settings which checks this
ETsamePr = False
if "ETsamePr" in binding:
if returnBool('ETsamePr'):
ETsamePr = True
if ETsamePr:
self.var.ETRef, MaskMapBoundary = readmeteodata(self.var.evaTMaps, dateVar['currDate'], addZeros=True, mapsscale=self.var.meteomapsscale)
self.var.ETRef = self.var.ETRef *self.var.DtDay * self.var.con_e
self.var.ETRef = self.downscaling2(self.var.ETRef, "downscale_wordclim_prec", self.var.wc2_prec, self.var.wc4_prec, downscale=0)
self.var.EWRef, MaskMapBoundary = readmeteodata(self.var.eva0Maps, dateVar['currDate'], addZeros=True, mapsscale=self.var.meteomapsscale)
self.var.EWRef = self.var.EWRef * self.var.DtDay * self.var.con_e
self.var.EWRef = self.downscaling2(self.var.EWRef, "downscale_wordclim_prec", self.var.wc2_prec, self.var.wc4_prec, downscale=0)
else:
self.var.ETRef, MaskMapBoundary = readmeteodata(self.var.evaTMaps, dateVar['currDate'], addZeros=True, mapsscale = True)
self.var.ETRef = self.var.ETRef *self.var.DtDay * self.var.con_e
self.var.EWRef, MaskMapBoundary = readmeteodata(self.var.eva0Maps, dateVar['currDate'], addZeros=True, mapsscale = True)
self.var.EWRef = self.var.EWRef * self.var.DtDay * self.var.con_e
# potential evaporation rate from water surface (conversion to [m] per time step)
# potential evaporation rate from a bare soil surface (conversion # to [m] per time step)
if Flags['calib']:
# if first clibration run, store all meteo data in a variable
if dateVar['curr'] == 1:
if self.var.includeGlaciers:
if not self.var.includeOnlyGlaciersMelt:
self.var.meteo = np.zeros([6, 1 + dateVar["intEnd"] - dateVar["intStart"], len(self.var.Precipitation)])
else:
self.var.meteo = np.zeros([5, 1 + dateVar["intEnd"] - dateVar["intStart"], len(self.var.Precipitation)])
else:
self.var.meteo = np.zeros([4,1 + dateVar["intEnd"] - dateVar["intStart"],len(self.var.Precipitation)])
self.var.meteo = np.zeros([4,1 + dateVar["intEnd"] - dateVar["intStart"],len(self.var.Precipitation)])
#self.var.ETRef_all,self.var.EWRef_all,self.var.Tavg_all, self.var.Precipitation_all = [],[],[],[]
no = dateVar['curr'] -1
self.var.meteo[0,no] = self.var.Precipitation
self.var.meteo[1,no] = self.var.Tavg
self.var.meteo[2,no] = self.var.ETRef
self.var.meteo[3,no] = self.var.EWRef
if self.var.includeGlaciers:
self.var.meteo[4, no] = self.var.GlacierMelt
if not self.var.includeOnlyGlaciersMelt:
self.var.meteo[5, no] = self.var.GlacierRain
#self.var.ETRef_all.append(self.var.ETRef)
#self.var.EWRef_all.append(self.var.EWRef)
#self.var.Tavg_all.append(self.var.Tavg)
#self.var.Precipitation_all.append(self.var.Precipitation)
ii =1