Index: wflow-py/wflow/pcrglobwb/virtualOS.py =================================================================== diff -u -r4143969ebd6b276284ea6e4198ee693bd6ab506a -rea74d7d019652362ccf1714f1c6448955e24b49d --- wflow-py/wflow/pcrglobwb/virtualOS.py (.../virtualOS.py) (revision 4143969ebd6b276284ea6e4198ee693bd6ab506a) +++ wflow-py/wflow/pcrglobwb/virtualOS.py (.../virtualOS.py) (revision ea74d7d019652362ccf1714f1c6448955e24b49d) @@ -79,9 +79,7 @@ + str(ncFile) ) - logger.debug('Check whether the variable: '+str(varName)+' is defined in the file: '+str(ncFile)) - - if ncFile in list(filecache.keys()): + if ncFile in filecache.keys(): f = filecache[ncFile] # ~ print "Cached: ", ncFile else: @@ -90,10 +88,10 @@ # ~ print "New: ", ncFile varName = str(varName) - - return varName in list(f.variables.keys()) + return varName in f.variables.keys() + def netcdf2PCRobjCloneWithoutTime( ncFile, varName, @@ -114,7 +112,7 @@ # Only works if cells are 'square'. # Only works if cellsizeClone <= cellsizeInput # Get netCDF file and variable name: - if ncFile in list(filecache.keys()): + if ncFile in filecache.keys(): f = filecache[ncFile] # ~ print "Cached: ", ncFile else: @@ -221,12 +219,12 @@ # Only works if cells are 'square'. # Only works if cellsizeClone <= cellsizeInput # Get netCDF file and variable name: - - #~ print ncFile - - logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile)) - - if ncFile in list(filecache.keys()): + + # ~ print ncFile + + logger.debug("reading variable: " + str(varName) + " from the file: " + str(ncFile)) + + if ncFile in filecache.keys(): f = filecache[ncFile] # ~ print "Cached: ", ncFile else: @@ -569,12 +567,12 @@ # Only works if cells are 'square'. # Only works if cellsizeClone <= cellsizeInput # Get netCDF file and variable name: - - #~ print ncFile - - logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile)) - - if ncFile in list(filecache.keys()): + + # ~ print ncFile + + logger.debug("reading variable: " + str(varName) + " from the file: " + str(ncFile)) + + if ncFile in filecache.keys(): f = filecache[ncFile] # ~ print "Cached: ", ncFile else: @@ -1421,26 +1419,50 @@ return False - if err !=None or cOut == []: - print("Something wrong with mattattr in virtualOS, maybe clone Map does not exist ? ") +def getMapAttributesALL(cloneMap, arcDegree=True): + cOut, err = subprocess.Popen( + str("mapattr -p %s " % (cloneMap)), + stdout=subprocess.PIPE, + stderr=open(os.devnull), + shell=True, + ).communicate() + + if err != None or cOut == []: + print "Something wrong with mattattr in virtualOS, maybe clone Map does not exist ? " sys.exit() cellsize = float(cOut.split()[7]) - if arcDegree == True: cellsize = round(cellsize * 360000.)/360000. - mapAttr = {'cellsize': float(cellsize) ,\ - 'rows' : float(cOut.split()[3]) ,\ - 'cols' : float(cOut.split()[5]) ,\ - 'xUL' : float(cOut.split()[17]),\ - 'yUL' : float(cOut.split()[19])} - co = None; cOut = None; err = None - del co; del cOut; del err - n = gc.collect() ; del gc.garbage[:] ; n = None ; del n - return mapAttr + if arcDegree == True: + cellsize = round(cellsize * 360000.) / 360000. + mapAttr = { + "cellsize": float(cellsize), + "rows": float(cOut.split()[3]), + "cols": float(cOut.split()[5]), + "xUL": float(cOut.split()[17]), + "yUL": float(cOut.split()[19]), + } + co = None + cOut = None + err = None + del co + del cOut + del err + n = gc.collect() + del gc.garbage[:] + n = None + del n + return mapAttr -def getMapAttributes(cloneMap,attribute,arcDegree=True): - cOut,err = subprocess.Popen(str('mapattr -p %s ' %(cloneMap)), stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - #print cOut - if err !=None or cOut == []: - print("Something wrong with mattattr in virtualOS, maybe clone Map does not exist ? ") + +def getMapAttributes(cloneMap, attribute, arcDegree=True): + cOut, err = subprocess.Popen( + str("mapattr -p %s " % (cloneMap)), + stdout=subprocess.PIPE, + stderr=open(os.devnull), + shell=True, + ).communicate() + # print cOut + if err != None or cOut == []: + print "Something wrong with mattattr in virtualOS, maybe clone Map does not exist ? " sys.exit() # print cOut.split() co = None @@ -1551,36 +1573,38 @@ def secondsPerDay(): return float(3600 * 24) - -def getValDivZero(x,y,y_lim=smallNumber,z_def= 0.): - #-returns the result of a division that possibly involves a zero - # denominator; in which case, a default value is substituted: - # x/y= z in case y > y_lim, - # x/y= z_def in case y <= y_lim, where y_lim -> 0. - # z_def is set to zero if not otherwise specified - return pcr.ifthenelse(y > y_lim,x/pcr.max(y_lim,y),z_def) -def getValFloatDivZero(x,y,y_lim,z_def= 0.): - #-returns the result of a division that possibly involves a zero - # denominator; in which case, a default value is substituted: - # x/y= z in case y > y_lim, - # x/y= z_def in case y <= y_lim, where y_lim -> 0. - # z_def is set to zero if not otherwise specified - if y > y_lim: - return x / max(y_lim,y) - else: - return z_def +def getValDivZero(x, y, y_lim=smallNumber, z_def=0.): + # -returns the result of a division that possibly involves a zero + # denominator; in which case, a default value is substituted: + # x/y= z in case y > y_lim, + # x/y= z_def in case y <= y_lim, where y_lim -> 0. + # z_def is set to zero if not otherwise specified + return pcr.ifthenelse(y > y_lim, x / pcr.max(y_lim, y), z_def) -def retrieveMapValue(pcrX,coordinates): - #-retrieves values from a map and returns an array conform the IDs stored in properties - nrRows= coordinates.shape[0] - x= np.ones((nrRows))* MV - tmpIDArray= pcr.pcr2numpy(pcrX,MV) - for iCnt in range(nrRows): - row,col= coordinates[iCnt,:] - if row != MV and col != MV: - x[iCnt]= tmpIDArray[row,col] + +def getValFloatDivZero(x, y, y_lim, z_def=0.): + # -returns the result of a division that possibly involves a zero + # denominator; in which case, a default value is substituted: + # x/y= z in case y > y_lim, + # x/y= z_def in case y <= y_lim, where y_lim -> 0. + # z_def is set to zero if not otherwise specified + if y > y_lim: + return x / max(y_lim, y) + else: + return z_def + + +def retrieveMapValue(pcrX, coordinates): + # -retrieves values from a map and returns an array conform the IDs stored in properties + nrRows = coordinates.shape[0] + x = np.ones((nrRows)) * MV + tmpIDArray = pcr.pcr2numpy(pcrX, MV) + for iCnt in xrange(nrRows): + row, col = coordinates[iCnt, :] + if row != MV and col != MV: + x[iCnt] = tmpIDArray[row, col] return x @@ -1589,16 +1613,16 @@ if x.ndim == 1: nrRows = 1 - tempIDArray= pcr.pcr2numpy(pcrX,MV) - #print tempIDArray - temporary= tempIDArray - nrRows= coord.shape[0] - for iCnt in range(nrRows): - row,col= coord[iCnt,:] - if row != MV and col != MV: - tempIDArray[row,col]= (x[iCnt]) - # print iCnt,row,col,x[iCnt] - pcrX= pcr.numpy2pcr(pcr.Scalar,tempIDArray,MV) + tempIDArray = pcr.pcr2numpy(pcrX, MV) + # print tempIDArray + temporary = tempIDArray + nrRows = coord.shape[0] + for iCnt in xrange(nrRows): + row, col = coord[iCnt, :] + if row != MV and col != MV: + tempIDArray[row, col] = x[iCnt] + # print iCnt,row,col,x[iCnt] + pcrX = pcr.numpy2pcr(pcr.Scalar, tempIDArray, MV) return pcrX @@ -1771,7 +1795,7 @@ # if abs(a) > 1e-5 or abs(b) > 1e-5: # if abs(a) > 1e-4 or abs(b) > 1e-4: if abs(a) > threshold or abs(b) > threshold: - print("WBError %s Min %f Max %f Mean %f" %(processName,a,b,c)) + print "WBError %s Min %f Max %f Mean %f" % (processName, a, b, c) # if abs(inflow + deltaS - outflow) > 1e-5: # print "Water balance Error for %s on %s: in = %f\tout=%f\tdeltaS=%f\tBalance=%f" \ # %(processName,dateStr,inflow,outflow,deltaS,inflow + deltaS - outflow) @@ -1880,14 +1904,27 @@ ) # water balance check - if debug_water_balance and not isinstance(zone_area,type(None)): - waterBalanceCheck([pcr.cover(pcr.areatotal(cellAbstraction, allocation_zones)/zone_area, 0.0)],\ - [pcr.cover(pcr.areatotal(cellAllocation , allocation_zones)/zone_area, 0.0)],\ - [pcr.scalar(0.0)],\ - [pcr.scalar(0.0)],\ - 'abstraction - allocation per zone/segment (with high precision) - loop (power number): ' + str(power_number) ,\ - True,\ - extra_info_for_water_balance_reporting, threshold = 1e-5) + if debug_water_balance and not isinstance(zone_area, types.NoneType): + waterBalanceCheck( + [ + pcr.cover( + pcr.areatotal(cellAbstraction, allocation_zones) / zone_area, + 0.0, + ) + ], + [ + pcr.cover( + pcr.areatotal(cellAllocation, allocation_zones) / zone_area, 0.0 + ) + ], + [pcr.scalar(0.0)], + [pcr.scalar(0.0)], + "abstraction - allocation per zone/segment (with high precision) - loop (power number): " + + str(power_number), + True, + extra_info_for_water_balance_reporting, + threshold=1e-5, + ) # actual water abstraction and allocation in this current loop (power number) cell_abstrac_for_every_power_number[str(power_number)] = cellAbstraction @@ -1905,15 +1942,26 @@ sumCellAllocation += cell_allocat_for_every_power_number[str(power_number)] # water balance check - if debug_water_balance and not isinstance(zone_area,type(None)): - waterBalanceCheck([pcr.cover(pcr.areatotal(sumCellAbstraction, allocation_zones)/zone_area, 0.0)],\ - [pcr.cover(pcr.areatotal(sumCellAllocation , allocation_zones)/zone_area, 0.0)],\ - [pcr.scalar(0.0)],\ - [pcr.scalar(0.0)],\ - 'abstraction - allocation per zone/segment (with high precision) - sum after loop' ,\ - True,\ - extra_info_for_water_balance_reporting, threshold = 1e-5) - + if debug_water_balance and not isinstance(zone_area, types.NoneType): + waterBalanceCheck( + [ + pcr.cover( + pcr.areatotal(sumCellAbstraction, allocation_zones) / zone_area, 0.0 + ) + ], + [ + pcr.cover( + pcr.areatotal(sumCellAllocation, allocation_zones) / zone_area, 0.0 + ) + ], + [pcr.scalar(0.0)], + [pcr.scalar(0.0)], + "abstraction - allocation per zone/segment (with high precision) - sum after loop", + True, + extra_info_for_water_balance_reporting, + threshold=1e-5, + ) + return sumCellAbstraction, sumCellAllocation @@ -1933,7 +1981,7 @@ # demand volume in each cell (unit: m3) cellVolDemand = pcr.max(0.0, water_demand_volume) - if not isinstance(landmask, type(None)): + if not isinstance(landmask, types.NoneType): cellVolDemand = pcr.ifthen(landmask, pcr.cover(cellVolDemand, 0.0)) if ignore_small_values: # ignore small values to avoid runding error cellVolDemand = pcr.rounddown(pcr.max(0.0, water_demand_volume)) @@ -1945,16 +1993,16 @@ # total available water volume in each cell cellAvlWater = pcr.max(0.0, available_water_volume) - if not isinstance(landmask, type(None)): + if not isinstance(landmask, types.NoneType): cellAvlWater = pcr.ifthen(landmask, pcr.cover(cellAvlWater, 0.0)) if ignore_small_values: # ignore small values to avoid runding error cellAvlWater = pcr.rounddown(pcr.max(0.00, available_water_volume)) else: cellAvlWater = pcr.max(0.0, available_water_volume) # total available water volume in each zone/segment (unit: m3) - # - to minimize numerical errors, separating cellAvlWater - if not isinstance(high_volume_treshold,type(None)): + # - to minimize numerical errors, separating cellAvlWater + if not isinstance(high_volume_treshold, types.NoneType): # mask: 0 for small volumes ; 1 for large volumes (e.g. in lakes and reservoirs) mask = pcr.cover( pcr.ifthen(cellAvlWater > high_volume_treshold, pcr.boolean(1)), @@ -1980,8 +2028,8 @@ cellAbstraction = pcr.min(cellAbstraction, cellAvlWater) if ignore_small_values: # ignore small values to avoid runding error cellAbstraction = pcr.rounddown(pcr.max(0.00, cellAbstraction)) - # to minimize numerical errors, separating cellAbstraction - if not isinstance(high_volume_treshold,type(None)): + # to minimize numerical errors, separating cellAbstraction + if not isinstance(high_volume_treshold, types.NoneType): # mask: 0 for small volumes ; 1 for large volumes (e.g. in lakes and reservoirs) mask = pcr.cover( pcr.ifthen(cellAbstraction > high_volume_treshold, pcr.boolean(1)), @@ -2021,28 +2069,40 @@ - pcr.areatotal(cellAllocation, allocation_zones), ) remainingCellDemand = pcr.max(0.0, cellVolDemand - cellAllocation) - cellAllocation += zoneDeficitAllocation * getValDivZero(\ - remainingCellDemand, - pcr.areatotal(remainingCellDemand, allocation_zones), - smallNumber) - - if debug_water_balance and not isinstance(zone_area,type(None)): + cellAllocation += zoneDeficitAllocation * getValDivZero( + remainingCellDemand, + pcr.areatotal(remainingCellDemand, allocation_zones), + smallNumber, + ) - waterBalanceCheck([pcr.cover(pcr.areatotal(cellAbstraction, allocation_zones)/zone_area, 0.0)],\ - [pcr.cover(pcr.areatotal(cellAllocation , allocation_zones)/zone_area, 0.0)],\ - [pcr.scalar(0.0)],\ - [pcr.scalar(0.0)],\ - 'abstraction - allocation per zone/segment (PS: Error here may be caused by rounding error.)' ,\ - True,\ - extra_info_for_water_balance_reporting,threshold=1e-4) - + if debug_water_balance and not isinstance(zone_area, types.NoneType): + + waterBalanceCheck( + [ + pcr.cover( + pcr.areatotal(cellAbstraction, allocation_zones) / zone_area, 0.0 + ) + ], + [ + pcr.cover( + pcr.areatotal(cellAllocation, allocation_zones) / zone_area, 0.0 + ) + ], + [pcr.scalar(0.0)], + [pcr.scalar(0.0)], + "abstraction - allocation per zone/segment (PS: Error here may be caused by rounding error.)", + True, + extra_info_for_water_balance_reporting, + threshold=1e-4, + ) + return cellAbstraction, cellAllocation def findLastYearInNCFile(ncFile): # open a netcdf file: - if ncFile in list(filecache.keys()): + if ncFile in filecache.keys(): f = filecache[ncFile] else: f = nc.Dataset(ncFile)