Index: wflow-py/wflow/pcrglobwb/virtualOS.py =================================================================== diff -u -r1738c121cfc83b76d915b8ed55747cb1bb9a2f64 -r3e39e84af48f1bcb5ec0d243748147be223674f2 --- wflow-py/wflow/pcrglobwb/virtualOS.py (.../virtualOS.py) (revision 1738c121cfc83b76d915b8ed55747cb1bb9a2f64) +++ wflow-py/wflow/pcrglobwb/virtualOS.py (.../virtualOS.py) (revision 3e39e84af48f1bcb5ec0d243748147be223674f2) @@ -46,1106 +46,1498 @@ from wflow.wflow_lib import getgridparams, getrows import logging -logger = logging.getLogger('wflow_pcrglobwb') +logger = logging.getLogger("wflow_pcrglobwb") -# file cache to minimize/reduce opening/closing files. + +# file cache to minimize/reduce opening/closing files. filecache = dict() # Global variables: MV = 1e20 smallNumber = 1E-39 # tuple of netcdf file suffixes (extensions) that can be used: -netcdf_suffixes = ('.nc4','.nc') +netcdf_suffixes = (".nc4", ".nc") + def getFileList(inputDir, filePattern): - '''creates a dictionary of files meeting the pattern specified''' - fileNameList = glob.glob(os.path.join(inputDir, filePattern)) - ll= {} - for fileName in fileNameList: - ll[os.path.split(fileName)[-1]]= fileName - return ll + """creates a dictionary of files meeting the pattern specified""" + fileNameList = glob.glob(os.path.join(inputDir, filePattern)) + ll = {} + for fileName in fileNameList: + ll[os.path.split(fileName)[-1]] = fileName + return ll -def checkVariableInNC(ncFile,varName): - logger.debug('Check whether the variable: '+str(varName)+' is defined in the file: '+str(ncFile)) - +def checkVariableInNC(ncFile, varName): + + logger.debug( + "Check whether the variable: " + + str(varName) + + " is defined in the file: " + + str(ncFile) + ) + if ncFile in list(filecache.keys()): f = filecache[ncFile] - #~ print "Cached: ", ncFile + # ~ print "Cached: ", ncFile else: f = nc.Dataset(ncFile) filecache[ncFile] = f - #~ print "New: ", ncFile - + # ~ print "New: ", ncFile + varName = str(varName) - + return varName in list(f.variables.keys()) -def netcdf2PCRobjCloneWithoutTime(ncFile, varName, - cloneMapFileName = None,\ - LatitudeLongitude = True,\ - specificFillValue = None,\ - absolutePath = None): - - if absolutePath != None: ncFile = getFullPath(ncFile, absolutePath) - - logger.debug('reading variable: '+str(varName)+' from the file: '+str(ncFile)) - - # +def netcdf2PCRobjCloneWithoutTime( + ncFile, + varName, + cloneMapFileName=None, + LatitudeLongitude=True, + specificFillValue=None, + absolutePath=None, +): + + if absolutePath != None: + ncFile = getFullPath(ncFile, absolutePath) + + logger.debug("reading variable: " + str(varName) + " from the file: " + str(ncFile)) + + # # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file. # --- with clone checking # Only works if cells are 'square'. # Only works if cellsizeClone <= cellsizeInput # Get netCDF file and variable name: if ncFile in list(filecache.keys()): f = filecache[ncFile] - #~ print "Cached: ", ncFile + # ~ print "Cached: ", ncFile else: f = nc.Dataset(ncFile) filecache[ncFile] = f - #~ print "New: ", ncFile - - #print ncFile - #f = nc.Dataset(ncFile) + # ~ print "New: ", ncFile + + # print ncFile + # f = nc.Dataset(ncFile) varName = str(varName) - + if LatitudeLongitude == True: try: - f.variables['lat'] = f.variables['latitude'] - f.variables['lon'] = f.variables['longitude'] + f.variables["lat"] = f.variables["latitude"] + f.variables["lon"] = f.variables["longitude"] except: pass - -# sameClone = True -# # check whether clone and input maps have the same attributes: -# if cloneMapFileName != None: -# # get the attributes of cloneMap -# #attributeClone = getMapAttributesALL(cloneMapFileName) -# #cellsizeClone = attributeClone['cellsize'] -# #rowsClone = attributeClone['rows'] -# #colsClone = attributeClone['cols'] -# #xULClone = attributeClone['xUL'] -# #yULClone = attributeClone['yUL'] -# attributeClone = getgridparams() -# cellsizeClone = attributeClone[2] -# rowsClone = attributeClone[4] -# colsClone = attributeClone[5] -# xULClone = attributeClone[0] - 0.5*cellsizeClone -# yULClone = attributeClone[1] + 0.5*cellsizeClone -# # get the attributes of input (netCDF) -# cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] -# cellsizeInput = float(cellsizeInput) -# rowsInput = len(f.variables['lat']) -# colsInput = len(f.variables['lon']) -# xULInput = f.variables['lon'][0]-0.5*cellsizeInput -# yULInput = f.variables['lat'][0]+0.5*cellsizeInput -# # check whether both maps have the same attributes -# if cellsizeClone != cellsizeInput: sameClone = False -# if rowsClone != rowsInput: sameClone = False -# if colsClone != colsInput: sameClone = False -# if xULClone != xULInput: sameClone = False -# if yULClone != yULInput: sameClone = False -# - cropData = f.variables[varName][:,:] # still original data -# factor = 1 # needed in regridData2FinerGrid -# if sameClone == False: -# # crop to cloneMap: -# minX = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX) -# xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0]) -# xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) -# minY = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY) -# yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0]) -# yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) -# cropData = f.variables[varName][yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] -# factor = int(round(float(cellsizeInput)/float(cellsizeClone))) -# -# if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone))) - + + # sameClone = True + # # check whether clone and input maps have the same attributes: + # if cloneMapFileName != None: + # # get the attributes of cloneMap + # #attributeClone = getMapAttributesALL(cloneMapFileName) + # #cellsizeClone = attributeClone['cellsize'] + # #rowsClone = attributeClone['rows'] + # #colsClone = attributeClone['cols'] + # #xULClone = attributeClone['xUL'] + # #yULClone = attributeClone['yUL'] + # attributeClone = getgridparams() + # cellsizeClone = attributeClone[2] + # rowsClone = attributeClone[4] + # colsClone = attributeClone[5] + # xULClone = attributeClone[0] - 0.5*cellsizeClone + # yULClone = attributeClone[1] + 0.5*cellsizeClone + # # get the attributes of input (netCDF) + # cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] + # cellsizeInput = float(cellsizeInput) + # rowsInput = len(f.variables['lat']) + # colsInput = len(f.variables['lon']) + # xULInput = f.variables['lon'][0]-0.5*cellsizeInput + # yULInput = f.variables['lat'][0]+0.5*cellsizeInput + # # check whether both maps have the same attributes + # if cellsizeClone != cellsizeInput: sameClone = False + # if rowsClone != rowsInput: sameClone = False + # if colsClone != colsInput: sameClone = False + # if xULClone != xULInput: sameClone = False + # if yULClone != yULInput: sameClone = False + # + cropData = f.variables[varName][:, :] # still original data + # factor = 1 # needed in regridData2FinerGrid + # if sameClone == False: + # # crop to cloneMap: + # minX = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX) + # xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0]) + # xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) + # minY = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY) + # yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0]) + # yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) + # cropData = f.variables[varName][yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] + # factor = int(round(float(cellsizeInput)/float(cellsizeClone))) + # + # if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone))) + # convert to PCR object and close f if specificFillValue != None: - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ - #regridData2FinerGrid(factor,cropData,MV), \ - float(specificFillValue)) + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(specificFillValue), + ) else: - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ - #regridData2FinerGrid(factor,cropData,MV), \ - float(f.variables[varName]._FillValue)) - - #~ # debug: - #~ pcr.report(outPCR,"tmp.map") - #~ print(varName) - #~ os.system('aguila tmp.map') - - #f.close(); - f = None ; cropData = None + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(f.variables[varName]._FillValue), + ) + + # ~ # debug: + # ~ pcr.report(outPCR,"tmp.map") + # ~ print(varName) + # ~ os.system('aguila tmp.map') + + # f.close(); + f = None + cropData = None # PCRaster object - return (outPCR) + return outPCR -def netcdf2PCRobjClone(ncFile,varName,dateInput,\ - useDoy = None, - cloneMapFileName = None,\ - LatitudeLongitude = True,\ - specificFillValue = None): - # + +def netcdf2PCRobjClone( + ncFile, + varName, + dateInput, + useDoy=None, + cloneMapFileName=None, + LatitudeLongitude=True, + specificFillValue=None, +): + # # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file. # --- with clone checking # 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)) - + + # ~ print ncFile + + logger.debug("reading variable: " + str(varName) + " from the file: " + str(ncFile)) + if ncFile in list(filecache.keys()): f = filecache[ncFile] - #~ print "Cached: ", ncFile + # ~ print "Cached: ", ncFile else: f = nc.Dataset(ncFile) filecache[ncFile] = f - #~ print "New: ", ncFile - + # ~ print "New: ", ncFile + varName = str(varName) - + if LatitudeLongitude == True: try: - f.variables['lat'] = f.variables['latitude'] - f.variables['lon'] = f.variables['longitude'] + f.variables["lat"] = f.variables["latitude"] + f.variables["lon"] = f.variables["longitude"] except: pass - - if varName == "evapotranspiration": + + if varName == "evapotranspiration": try: - f.variables['evapotranspiration'] = f.variables['referencePotET'] + f.variables["evapotranspiration"] = f.variables["referencePotET"] except: pass - if varName == "kc": # the variable name in PCR-GLOBWB - try: - f.variables['kc'] = \ - f.variables['Cropcoefficient'] # the variable name in the netcdf file - except: - pass + if varName == "kc": # the variable name in PCR-GLOBWB + try: + f.variables["kc"] = f.variables[ + "Cropcoefficient" + ] # the variable name in the netcdf file + except: + pass - if varName == "interceptCapInput": # the variable name in PCR-GLOBWB - try: - f.variables['interceptCapInput'] = \ - f.variables['Interceptioncapacity'] # the variable name in the netcdf file - except: - pass + if varName == "interceptCapInput": # the variable name in PCR-GLOBWB + try: + f.variables["interceptCapInput"] = f.variables[ + "Interceptioncapacity" + ] # the variable name in the netcdf file + except: + pass - if varName == "coverFractionInput": # the variable name in PCR-GLOBWB - try: - f.variables['coverFractionInput'] = \ - f.variables['Coverfraction'] # the variable name in the netcdf file - except: - pass + if varName == "coverFractionInput": # the variable name in PCR-GLOBWB + try: + f.variables["coverFractionInput"] = f.variables[ + "Coverfraction" + ] # the variable name in the netcdf file + except: + pass - if varName == "fracVegCover": # the variable name in PCR-GLOBWB - try: - f.variables['fracVegCover'] = \ - f.variables['vegetation_fraction'] # the variable name in the netcdf file - except: - pass + if varName == "fracVegCover": # the variable name in PCR-GLOBWB + try: + f.variables["fracVegCover"] = f.variables[ + "vegetation_fraction" + ] # the variable name in the netcdf file + except: + pass - if varName == "minSoilDepthFrac": # the variable name in PCR-GLOBWB - try: - f.variables['minSoilDepthFrac'] = \ - f.variables['minRootDepthFraction'] # the variable name in the netcdf file - except: - pass + if varName == "minSoilDepthFrac": # the variable name in PCR-GLOBWB + try: + f.variables["minSoilDepthFrac"] = f.variables[ + "minRootDepthFraction" + ] # the variable name in the netcdf file + except: + pass - if varName == "maxSoilDepthFrac": # the variable name in PCR-GLOBWB - try: - f.variables['maxSoilDepthFrac'] = \ - f.variables['maxRootDepthFraction'] # the variable name in the netcdf file - except: - pass + if varName == "maxSoilDepthFrac": # the variable name in PCR-GLOBWB + try: + f.variables["maxSoilDepthFrac"] = f.variables[ + "maxRootDepthFraction" + ] # the variable name in the netcdf file + except: + pass - if varName == "arnoBeta": # the variable name in PCR-GLOBWB - try: - f.variables['arnoBeta'] = \ - f.variables['arnoSchemeBeta'] # the variable name in the netcdf file - except: - pass + if varName == "arnoBeta": # the variable name in PCR-GLOBWB + try: + f.variables["arnoBeta"] = f.variables[ + "arnoSchemeBeta" + ] # the variable name in the netcdf file + except: + pass # date date = dateInput - if useDoy == "Yes": - logger.debug('Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)') + if useDoy == "Yes": + logger.debug( + "Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)" + ) idx = int(dateInput) - 1 - elif useDoy == "month": # PS: WE NEED THIS ONE FOR NETCDF FILES that contain only 12 monthly values (e.g. cropCoefficientWaterNC). - logger.debug('Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)') + elif ( + useDoy == "month" + ): # PS: WE NEED THIS ONE FOR NETCDF FILES that contain only 12 monthly values (e.g. cropCoefficientWaterNC). + logger.debug( + "Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)" + ) # make sure that date is in the correct format - if isinstance(date, str) == True: date = \ - datetime.datetime.strptime(str(date),'%Y-%m-%d') + if isinstance(date, str) == True: + date = datetime.datetime.strptime(str(date), "%Y-%m-%d") idx = int(date.month) - 1 else: # make sure that date is in the correct format - if isinstance(date, str) == True: date = \ - datetime.datetime.strptime(str(date),'%Y-%m-%d') - date = datetime.datetime(date.year,date.month,date.day) + if isinstance(date, str) == True: + date = datetime.datetime.strptime(str(date), "%Y-%m-%d") + date = datetime.datetime(date.year, date.month, date.day) if useDoy == "yearly": - date = datetime.datetime(date.year,int(1),int(1)) + date = datetime.datetime(date.year, int(1), int(1)) if useDoy == "monthly": - date = datetime.datetime(date.year,date.month,int(1)) + date = datetime.datetime(date.year, date.month, int(1)) if useDoy == "yearly" or useDoy == "monthly" or useDoy == "daily_seasonal": # if the desired year is not available, use the first year or the last year that is available - first_year_in_nc_file = findFirstYearInNCTime(f.variables['time']) - last_year_in_nc_file = findLastYearInNCTime(f.variables['time']) + first_year_in_nc_file = findFirstYearInNCTime(f.variables["time"]) + last_year_in_nc_file = findLastYearInNCTime(f.variables["time"]) # - if date.year < first_year_in_nc_file: - if date.day == 29 and date.month == 2 and calendar.isleap(date.year) and calendar.isleap(first_year_in_nc_file) == False: + if date.year < first_year_in_nc_file: + if ( + date.day == 29 + and date.month == 2 + and calendar.isleap(date.year) + and calendar.isleap(first_year_in_nc_file) == False + ): date = datetime.datetime(first_year_in_nc_file, date.month, 28) else: - date = datetime.datetime(first_year_in_nc_file, date.month, date.day) - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(dateInput)+" is NOT available. " - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used." - # msg += "\n" + date = datetime.datetime( + first_year_in_nc_file, date.month, date.day + ) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += "The date " + str(dateInput) + " is NOT available. " + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is used." + ) + # msg += "\n" logger.warning(msg) - if date.year > last_year_in_nc_file: - if date.day == 29 and date.month == 2 and calendar.isleap(date.year) and calendar.isleap(last_year_in_nc_file) == False: + if date.year > last_year_in_nc_file: + if ( + date.day == 29 + and date.month == 2 + and calendar.isleap(date.year) + and calendar.isleap(last_year_in_nc_file) == False + ): date = datetime.datetime(last_year_in_nc_file, date.month, 28) else: date = datetime.datetime(last_year_in_nc_file, date.month, date.day) - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(dateInput)+" is NOT available. " - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used." - #msg += "\n" + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += "The date " + str(dateInput) + " is NOT available. " + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is used." + ) + # msg += "\n" logger.warning(msg) try: - idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \ - select ='exact') - msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is available. The 'exact' option is used while selecting netcdf time." + idx = nc.date2index( + date, + f.variables["time"], + calendar=f.variables["time"].calendar, + select="exact", + ) + msg = ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is available. The 'exact' option is used while selecting netcdf time." + ) logger.debug(msg) except: - msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'exact' option CANNOT be used while selecting netcdf time." + msg = ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is NOT available. The 'exact' option CANNOT be used while selecting netcdf time." + ) logger.debug(msg) - try: - idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \ - select = 'before') - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'before' option is used while selecting netcdf time." - #msg += "\n" + try: + idx = nc.date2index( + date, + f.variables["time"], + calendar=f.variables["time"].calendar, + select="before", + ) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is NOT available. The 'before' option is used while selecting netcdf time." + ) + # msg += "\n" except: - idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \ - select = 'after') - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'after' option is used while selecting netcdf time." - #msg += "\n" + idx = nc.date2index( + date, + f.variables["time"], + calendar=f.variables["time"].calendar, + select="after", + ) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is NOT available. The 'after' option is used while selecting netcdf time." + ) + # msg += "\n" logger.warning(msg) - - idx = int(idx) - logger.debug('Using the date index '+str(idx)) -# sameClone = True -# # check whether clone and input maps have the same attributes: -# if cloneMapFileName != None: -# # get the attributes of cloneMap -# #attributeClone = getMapAttributesALL(cloneMapFileName) -# #cellsizeClone = attributeClone['cellsize'] -# #rowsClone = attributeClone['rows'] -# #colsClone = attributeClone['cols'] -# #xULClone = attributeClone['xUL'] -# #yULClone = attributeClone['yUL'] -# attributeClone = getgridparams() -# cellsizeClone = attributeClone[2] -# rowsClone = attributeClone[4] -# colsClone = attributeClone[5] -# xULClone = attributeClone[0] - 0.5*cellsizeClone -# yULClone = attributeClone[1] + 0.5*cellsizeClone -# # get the attributes of input (netCDF) -# cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] -# cellsizeInput = float(cellsizeInput) -# rowsInput = len(f.variables['lat']) -# colsInput = len(f.variables['lon']) -# xULInput = f.variables['lon'][0]-0.5*cellsizeInput -# yULInput = f.variables['lat'][0]+0.5*cellsizeInput -# # check whether both maps have the same attributes -# if cellsizeClone != cellsizeInput: sameClone = False -# if rowsClone != rowsInput: sameClone = False -# if colsClone != colsInput: sameClone = False -# if xULClone != xULInput: sameClone = False -# if yULClone != yULInput: sameClone = False + idx = int(idx) + logger.debug("Using the date index " + str(idx)) - cropData = f.variables[varName][int(idx),:,:] # still original data -# factor = 1 # needed in regridData2FinerGrid + # sameClone = True + # # check whether clone and input maps have the same attributes: + # if cloneMapFileName != None: + # # get the attributes of cloneMap + # #attributeClone = getMapAttributesALL(cloneMapFileName) + # #cellsizeClone = attributeClone['cellsize'] + # #rowsClone = attributeClone['rows'] + # #colsClone = attributeClone['cols'] + # #xULClone = attributeClone['xUL'] + # #yULClone = attributeClone['yUL'] + # attributeClone = getgridparams() + # cellsizeClone = attributeClone[2] + # rowsClone = attributeClone[4] + # colsClone = attributeClone[5] + # xULClone = attributeClone[0] - 0.5*cellsizeClone + # yULClone = attributeClone[1] + 0.5*cellsizeClone + # # get the attributes of input (netCDF) + # cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] + # cellsizeInput = float(cellsizeInput) + # rowsInput = len(f.variables['lat']) + # colsInput = len(f.variables['lon']) + # xULInput = f.variables['lon'][0]-0.5*cellsizeInput + # yULInput = f.variables['lat'][0]+0.5*cellsizeInput + # # check whether both maps have the same attributes + # if cellsizeClone != cellsizeInput: sameClone = False + # if rowsClone != rowsInput: sameClone = False + # if colsClone != colsInput: sameClone = False + # if xULClone != xULInput: sameClone = False + # if yULClone != yULInput: sameClone = False -# if sameClone == False: -# -# logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone)) -# # crop to cloneMap: -# #~ xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0]) -# minX = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX) -# xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0]) -# xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) -# #~ yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0]) -# minY = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY) -# yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0]) -# yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) -# cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] -# -# factor = int(round(float(cellsizeInput)/float(cellsizeClone))) -# if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone))) + cropData = f.variables[varName][int(idx), :, :] # still original data + # factor = 1 # needed in regridData2FinerGrid + # if sameClone == False: + # + # logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone)) + # # crop to cloneMap: + # #~ xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0]) + # minX = min(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX) + # xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0]) + # xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) + # #~ yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0]) + # minY = min(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY) + # yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0]) + # yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) + # cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] + # + # factor = int(round(float(cellsizeInput)/float(cellsizeClone))) + # if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone))) + # convert to PCR object and close f if specificFillValue != None: - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ -# regridData2FinerGrid(factor,cropData,MV), \ - float(specificFillValue)) + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(specificFillValue), + ) else: - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ -# regridData2FinerGrid(factor,cropData,MV), \ - float(f.variables[varName]._FillValue)) - - #f.close(); - f = None ; cropData = None + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(f.variables[varName]._FillValue), + ) + + # f.close(); + f = None + cropData = None # PCRaster object - return (outPCR) + return outPCR -def netcdf2PCRobjCloneJOYCE(ncFile,varName,dateInput,\ - useDoy = None, - cloneMapFileName = None,\ - LatitudeLongitude = True,\ - specificFillValue = None): - # + +def netcdf2PCRobjCloneJOYCE( + ncFile, + varName, + dateInput, + useDoy=None, + cloneMapFileName=None, + LatitudeLongitude=True, + specificFillValue=None, +): + # # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file. # --- with clone checking # 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)) - + + # ~ print ncFile + + logger.debug("reading variable: " + str(varName) + " from the file: " + str(ncFile)) + if ncFile in list(filecache.keys()): f = filecache[ncFile] - #~ print "Cached: ", ncFile + # ~ print "Cached: ", ncFile else: f = nc.Dataset(ncFile) filecache[ncFile] = f - #~ print "New: ", ncFile - + # ~ print "New: ", ncFile + varName = str(varName) - + if LatitudeLongitude == True: try: - f.variables['lat'] = f.variables['latitude'] - f.variables['lon'] = f.variables['longitude'] + f.variables["lat"] = f.variables["latitude"] + f.variables["lon"] = f.variables["longitude"] except: pass - - if varName == "evapotranspiration": + + if varName == "evapotranspiration": try: - f.variables['evapotranspiration'] = f.variables['referencePotET'] + f.variables["evapotranspiration"] = f.variables["referencePotET"] except: pass - if varName == "kc": # the variable name in PCR-GLOBWB - try: - f.variables['kc'] = \ - f.variables['Cropcoefficient'] # the variable name in the netcdf file - except: - pass + if varName == "kc": # the variable name in PCR-GLOBWB + try: + f.variables["kc"] = f.variables[ + "Cropcoefficient" + ] # the variable name in the netcdf file + except: + pass - if varName == "interceptCapInput": # the variable name in PCR-GLOBWB - try: - f.variables['interceptCapInput'] = \ - f.variables['Interceptioncapacity'] # the variable name in the netcdf file - except: - pass + if varName == "interceptCapInput": # the variable name in PCR-GLOBWB + try: + f.variables["interceptCapInput"] = f.variables[ + "Interceptioncapacity" + ] # the variable name in the netcdf file + except: + pass - if varName == "coverFractionInput": # the variable name in PCR-GLOBWB - try: - f.variables['coverFractionInput'] = \ - f.variables['Coverfraction'] # the variable name in the netcdf file - except: - pass + if varName == "coverFractionInput": # the variable name in PCR-GLOBWB + try: + f.variables["coverFractionInput"] = f.variables[ + "Coverfraction" + ] # the variable name in the netcdf file + except: + pass - if varName == "fracVegCover": # the variable name in PCR-GLOBWB - try: - f.variables['fracVegCover'] = \ - f.variables['vegetation_fraction'] # the variable name in the netcdf file - except: - pass + if varName == "fracVegCover": # the variable name in PCR-GLOBWB + try: + f.variables["fracVegCover"] = f.variables[ + "vegetation_fraction" + ] # the variable name in the netcdf file + except: + pass - if varName == "arnoBeta": # the variable name in PCR-GLOBWB - try: - f.variables['arnoBeta'] = \ - f.variables['arnoSchemeBeta'] # the variable name in the netcdf file - except: - pass + if varName == "arnoBeta": # the variable name in PCR-GLOBWB + try: + f.variables["arnoBeta"] = f.variables[ + "arnoSchemeBeta" + ] # the variable name in the netcdf file + except: + pass # date date = dateInput - if useDoy == "Yes": - logger.debug('Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)') + if useDoy == "Yes": + logger.debug( + "Finding the date based on the given climatology doy index (1 to 366, or index 0 to 365)" + ) idx = int(dateInput) - 1 else: # make sure that date is in the correct format - if isinstance(date, str) == True: date = \ - datetime.datetime.strptime(str(date),'%Y-%m-%d') - date = datetime.datetime(date.year,date.month,date.day) + if isinstance(date, str) == True: + date = datetime.datetime.strptime(str(date), "%Y-%m-%d") + date = datetime.datetime(date.year, date.month, date.day) if useDoy == "month": - logger.debug('Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)') + logger.debug( + "Finding the date based on the given climatology month index (1 to 12, or index 0 to 11)" + ) idx = int(date.month) - 1 if useDoy == "yearly": - date = datetime.datetime(date.year,int(1),int(1)) + date = datetime.datetime(date.year, int(1), int(1)) if useDoy == "monthly": - date = datetime.datetime(date.year,date.month,int(1)) + date = datetime.datetime(date.year, date.month, int(1)) if useDoy == "yearly" or useDoy == "monthly" or useDoy == "daily_seasonal": # if the desired year is not available, use the first year or the last year that is available - first_year_in_nc_file = findFirstYearInNCTime(f.variables['time']) - last_year_in_nc_file = findLastYearInNCTime(f.variables['time']) + first_year_in_nc_file = findFirstYearInNCTime(f.variables["time"]) + last_year_in_nc_file = findLastYearInNCTime(f.variables["time"]) # - if date.year < first_year_in_nc_file: - date = datetime.datetime(first_year_in_nc_file,date.month,date.day) - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(dateInput)+" is NOT available. " - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used." - #msg += "\n" + if date.year < first_year_in_nc_file: + date = datetime.datetime(first_year_in_nc_file, date.month, date.day) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += "The date " + str(dateInput) + " is NOT available. " + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is used." + ) + # msg += "\n" logger.warning(msg) - if date.year > last_year_in_nc_file: - date = datetime.datetime(last_year_in_nc_file,date.month,date.day) - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(dateInput)+" is NOT available. " - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is used." - #msg += "\n" + if date.year > last_year_in_nc_file: + date = datetime.datetime(last_year_in_nc_file, date.month, date.day) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += "The date " + str(dateInput) + " is NOT available. " + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is used." + ) + # msg += "\n" logger.warning(msg) try: - idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \ - select ='exact') - msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is available. The 'exact' option is used while selecting netcdf time." + idx = nc.date2index( + date, + f.variables["time"], + calendar=f.variables["time"].calendar, + select="exact", + ) + msg = ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is available. The 'exact' option is used while selecting netcdf time." + ) logger.debug(msg) except: - msg = "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'exact' option CANNOT be used while selecting netcdf time." + msg = ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is NOT available. The 'exact' option CANNOT be used while selecting netcdf time." + ) logger.debug(msg) - try: - idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \ - select = 'before') - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'before' option is used while selecting netcdf time." - #msg += "\n" + try: + idx = nc.date2index( + date, + f.variables["time"], + calendar=f.variables["time"].calendar, + select="before", + ) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is NOT available. The 'before' option is used while selecting netcdf time." + ) + # msg += "\n" except: - idx = nc.date2index(date, f.variables['time'], calendar = f.variables['time'].calendar, \ - select = 'after') - #msg = "\n" - msg = "WARNING related to the netcdf file: "+str(ncFile)+" ; variable: "+str(varName)+" !!!!!! "#+"\n" - msg += "The date "+str(date.year)+"-"+str(date.month)+"-"+str(date.day)+" is NOT available. The 'after' option is used while selecting netcdf time." - #msg += "\n" + idx = nc.date2index( + date, + f.variables["time"], + calendar=f.variables["time"].calendar, + select="after", + ) + # msg = "\n" + msg = ( + "WARNING related to the netcdf file: " + + str(ncFile) + + " ; variable: " + + str(varName) + + " !!!!!! " + ) # +"\n" + msg += ( + "The date " + + str(date.year) + + "-" + + str(date.month) + + "-" + + str(date.day) + + " is NOT available. The 'after' option is used while selecting netcdf time." + ) + # msg += "\n" logger.warning(msg) - - idx = int(idx) - logger.debug('Using the date index '+str(idx)) - cropData = f.variables[varName][int(idx),:,:].copy() # still original data -# factor = 1 # needed in regridData2FinerGrid - + idx = int(idx) + logger.debug("Using the date index " + str(idx)) + + cropData = f.variables[varName][int(idx), :, :].copy() # still original data + # factor = 1 # needed in regridData2FinerGrid + # store latitudes and longitudes to a new variable - latitude = f.variables['lat'] - longitude = f.variables['lon'] - + latitude = f.variables["lat"] + longitude = f.variables["lon"] + # check the orientation of the latitude and flip it if necessary we_have_to_flip = False - if (latitude[0]- latitude[1]) < 0.0: + if (latitude[0] - latitude[1]) < 0.0: we_have_to_flip = True latitude = np.flipud(latitude) - -# sameClone = True -# # check whether clone and input maps have the same attributes: -# if cloneMapFileName != None: -# # get the attributes of cloneMap -# #attributeClone = getMapAttributesALL(cloneMapFileName) -# #cellsizeClone = attributeClone['cellsize'] -# #rowsClone = attributeClone['rows'] -# #colsClone = attributeClone['cols'] -# #xULClone = attributeClone['xUL'] -# #yULClone = attributeClone['yUL'] -# attributeClone = getgridparams() -# cellsizeClone = attributeClone[2] -# rowsClone = attributeClone[4] -# colsClone = attributeClone[5] -# xULClone = attributeClone[0] - 0.5*cellsizeClone -# yULClone = attributeClone[1] + 0.5*cellsizeClone -# # get the attributes of input (netCDF) -# cellsizeInput = latitude[0]- latitude[1] -# cellsizeInput = float(cellsizeInput) -# rowsInput = len(latitude) -# colsInput = len(longitude) -# xULInput = longitude[0]-0.5*cellsizeInput -# yULInput = latitude[0] +0.5*cellsizeInput -# # check whether both maps have the same attributes -# if cellsizeClone != cellsizeInput: sameClone = False -# if rowsClone != rowsInput: sameClone = False -# if colsClone != colsInput: sameClone = False -# if xULClone != xULInput: sameClone = False -# if yULClone != yULInput: sameClone = False - # flip cropData if necessary - if we_have_to_flip: - #~ cropData = cropData[::-1,:] - #~ cropData = cropData[::-1,:].copy() + # sameClone = True + # # check whether clone and input maps have the same attributes: + # if cloneMapFileName != None: + # # get the attributes of cloneMap + # #attributeClone = getMapAttributesALL(cloneMapFileName) + # #cellsizeClone = attributeClone['cellsize'] + # #rowsClone = attributeClone['rows'] + # #colsClone = attributeClone['cols'] + # #xULClone = attributeClone['xUL'] + # #yULClone = attributeClone['yUL'] + # attributeClone = getgridparams() + # cellsizeClone = attributeClone[2] + # rowsClone = attributeClone[4] + # colsClone = attributeClone[5] + # xULClone = attributeClone[0] - 0.5*cellsizeClone + # yULClone = attributeClone[1] + 0.5*cellsizeClone + # # get the attributes of input (netCDF) + # cellsizeInput = latitude[0]- latitude[1] + # cellsizeInput = float(cellsizeInput) + # rowsInput = len(latitude) + # colsInput = len(longitude) + # xULInput = longitude[0]-0.5*cellsizeInput + # yULInput = latitude[0] +0.5*cellsizeInput + # # check whether both maps have the same attributes + # if cellsizeClone != cellsizeInput: sameClone = False + # if rowsClone != rowsInput: sameClone = False + # if colsClone != colsInput: sameClone = False + # if xULClone != xULInput: sameClone = False + # if yULClone != yULInput: sameClone = False - #~ cropData = np.flipud(cropData) + # flip cropData if necessary + if we_have_to_flip: + # ~ cropData = cropData[::-1,:] + # ~ cropData = cropData[::-1,:].copy() - #~ cropData = np.flipud(cropData) - #~ cropData = np.flipud(cropData).copy() + # ~ cropData = np.flipud(cropData) - #~ original = cropData.copy() -#~ - #~ print id(cropData) - #~ print id(original) + # ~ cropData = np.flipud(cropData) + # ~ cropData = np.flipud(cropData).copy() - #~ cropData = None - #~ del cropData - #~ cropData = np.flipud(original).copy() - - #~ print type(cropData) - - #~ cropData2 = cropData[::-1,:] - - #~ cropData = None - #~ cropData = original[::-1,:] - #~ cropData = cropData[::-1,:] + # ~ original = cropData.copy() + # ~ + # ~ print id(cropData) + # ~ print id(original) - cropData = cropData[::-1,:] - - #print type(cropData) + # ~ cropData = None + # ~ del cropData + # ~ cropData = np.flipud(original).copy() - #print "Test test tet" - #print id(cropData) - #~ print id(original) + # ~ print type(cropData) - #~ cropData = cropData[::-1,:].copy() + # ~ cropData2 = cropData[::-1,:] + # ~ cropData = None + # ~ cropData = original[::-1,:] + # ~ cropData = cropData[::-1,:] + + cropData = cropData[::-1, :] + + # print type(cropData) + + # print "Test test tet" + # print id(cropData) + # ~ print id(original) + + # ~ cropData = cropData[::-1,:].copy() + pcr_map = pcr.numpy2pcr(pcr.Scalar, cropData, -999.9) pcr.report(pcr_map, "test2.map") os.system("aguila test2.map") - -# if sameClone == False: -# -# logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone)) -# # crop to cloneMap: -# minX = min(abs(longitude[:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX) -# xIdxSta = int(np.where(abs(longitude[:] - (xULClone + 0.5*cellsizeInput)) == minX)[0]) -# xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) -# minY = min(abs(latitude[:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY) -# yIdxSta = int(np.where(abs(latitude[:] - (yULClone - 0.5*cellsizeInput)) == minY)[0]) -# yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) -# cropData = cropData[yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] -# -# factor = int(round(float(cellsizeInput)/float(cellsizeClone))) -# if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone))) + # if sameClone == False: + # + # logger.debug('Crop to the clone map with lower left corner (x,y): '+str(xULClone)+' , '+str(yULClone)) + # # crop to cloneMap: + # minX = min(abs(longitude[:] - (xULClone + 0.5*cellsizeInput))) # ; print(minX) + # xIdxSta = int(np.where(abs(longitude[:] - (xULClone + 0.5*cellsizeInput)) == minX)[0]) + # xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) + # minY = min(abs(latitude[:] - (yULClone - 0.5*cellsizeInput))) # ; print(minY) + # yIdxSta = int(np.where(abs(latitude[:] - (yULClone - 0.5*cellsizeInput)) == minY)[0]) + # yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) + # cropData = cropData[yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] + # + # factor = int(round(float(cellsizeInput)/float(cellsizeClone))) + # if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone))) + # convert to PCR object and close f if specificFillValue != None: - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ -# regridData2FinerGrid(factor,cropData,MV), \ - float(specificFillValue)) + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(specificFillValue), + ) else: - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ -# regridData2FinerGrid(factor,cropData,MV), \ - float(f.variables[varName]._FillValue)) - - #f.close(); - f = None ; cropData = None + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(f.variables[varName]._FillValue), + ) + + # f.close(); + f = None + cropData = None # PCRaster object - return (outPCR) + return outPCR -def netcdf2PCRobjCloneWindDist(ncFile,varName,dateInput,useDoy = None, - cloneMapFileName=None): +def netcdf2PCRobjCloneWindDist( + ncFile, varName, dateInput, useDoy=None, cloneMapFileName=None +): # EHS (02 SEP 2013): This is a special function made by Niko Wanders (for his DA framework). # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file. # --- with clone checking # Only works if cells are 'square'. # Only works if cellsizeClone <= cellsizeInput - + # Get netCDF file and variable name: f = nc.Dataset(ncFile) varName = str(varName) # date date = dateInput - if useDoy == "Yes": + if useDoy == "Yes": idx = dateInput - 1 else: - if isinstance(date, str) == True: date = \ - datetime.datetime.strptime(str(date),'%Y-%m-%d') - date = datetime.datetime(date.year,date.month,date.day) + if isinstance(date, str) == True: + date = datetime.datetime.strptime(str(date), "%Y-%m-%d") + date = datetime.datetime(date.year, date.month, date.day) # time index (in the netCDF file) - nctime = f.variables['time'] # A netCDF time variable object. - idx = nc.date2index(date, nctime, calendar=nctime.calendar, \ - select='exact') - idx = int(idx) + nctime = f.variables["time"] # A netCDF time variable object. + idx = nc.date2index(date, nctime, calendar=nctime.calendar, select="exact") + idx = int(idx) -# sameClone = True -# # check whether clone and input maps have the same attributes: -# if cloneMapFileName != None: -# # get the attributes of cloneMap -# #attributeClone = getMapAttributesALL(cloneMapFileName) -# #cellsizeClone = attributeClone['cellsize'] -# #rowsClone = attributeClone['rows'] -# #colsClone = attributeClone['cols'] -# #xULClone = attributeClone['xUL'] -# #yULClone = attributeClone['yUL'] -# attributeClone = getgridparams() -# cellsizeClone = attributeClone[2] -# rowsClone = attributeClone[4] -# colsClone = attributeClone[5] -# xULClone = attributeClone[0] - 0.5*cellsizeClone -# yULClone = attributeClone[1] + 0.5*cellsizeClone -# # get the attributes of input (netCDF) -# cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] -# cellsizeInput = float(cellsizeInput) -# rowsInput = len(f.variables['lat']) -# colsInput = len(f.variables['lon']) -# xULInput = f.variables['lon'][0]-0.5*cellsizeInput -# yULInput = f.variables['lat'][0]+0.5*cellsizeInput -# # check whether both maps have the same attributes -# if cellsizeClone != cellsizeInput: sameClone = False -# if rowsClone != rowsInput: sameClone = False -# if colsClone != colsInput: sameClone = False -# if xULClone != xULInput: sameClone = False -# if yULClone != yULInput: sameClone = False + # sameClone = True + # # check whether clone and input maps have the same attributes: + # if cloneMapFileName != None: + # # get the attributes of cloneMap + # #attributeClone = getMapAttributesALL(cloneMapFileName) + # #cellsizeClone = attributeClone['cellsize'] + # #rowsClone = attributeClone['rows'] + # #colsClone = attributeClone['cols'] + # #xULClone = attributeClone['xUL'] + # #yULClone = attributeClone['yUL'] + # attributeClone = getgridparams() + # cellsizeClone = attributeClone[2] + # rowsClone = attributeClone[4] + # colsClone = attributeClone[5] + # xULClone = attributeClone[0] - 0.5*cellsizeClone + # yULClone = attributeClone[1] + 0.5*cellsizeClone + # # get the attributes of input (netCDF) + # cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] + # cellsizeInput = float(cellsizeInput) + # rowsInput = len(f.variables['lat']) + # colsInput = len(f.variables['lon']) + # xULInput = f.variables['lon'][0]-0.5*cellsizeInput + # yULInput = f.variables['lat'][0]+0.5*cellsizeInput + # # check whether both maps have the same attributes + # if cellsizeClone != cellsizeInput: sameClone = False + # if rowsClone != rowsInput: sameClone = False + # if colsClone != colsInput: sameClone = False + # if xULClone != xULInput: sameClone = False + # if yULClone != yULInput: sameClone = False - cropData = f.variables[varName][int(idx),:,:] # still original data -# factor = 1 # needed in regridData2FinerGrid -# if sameClone == False: -# # crop to cloneMap: -# xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0]) -# xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) -# yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0]) -# yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) -# cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] -# factor = int(float(cellsizeInput)/float(cellsizeClone)) - + cropData = f.variables[varName][int(idx), :, :] # still original data + # factor = 1 # needed in regridData2FinerGrid + # if sameClone == False: + # # crop to cloneMap: + # xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0]) + # xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) + # yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0]) + # yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) + # cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] + # factor = int(float(cellsizeInput)/float(cellsizeClone)) + # convert to PCR object and close f - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ -# regridData2FinerGrid(factor,cropData,MV), \ - float(0.0)) - f.close(); - f = None ; cropData = None + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(0.0), + ) + f.close() + f = None + cropData = None # PCRaster object - return (outPCR) - -def netcdf2PCRobjCloneWind(ncFile,varName,dateInput,useDoy = None, - cloneMapFileName=None): + return outPCR + + +def netcdf2PCRobjCloneWind( + ncFile, varName, dateInput, useDoy=None, cloneMapFileName=None +): # EHS (02 SEP 2013): This is a special function made by Niko Wanders (for his DA framework). # EHS (19 APR 2013): To convert netCDF (tss) file to PCR file. # --- with clone checking # Only works if cells are 'square'. # Only works if cellsizeClone <= cellsizeInput - + # Get netCDF file and variable name: f = nc.Dataset(ncFile) varName = str(varName) # date date = dateInput - if useDoy == "Yes": + if useDoy == "Yes": idx = dateInput - 1 else: - if isinstance(date, str) == True: date = \ - datetime.datetime.strptime(str(date),'%Y-%m-%d') - date = datetime.datetime(date.year,date.month,date.day, 0, 0) + if isinstance(date, str) == True: + date = datetime.datetime.strptime(str(date), "%Y-%m-%d") + date = datetime.datetime(date.year, date.month, date.day, 0, 0) # time index (in the netCDF file) - nctime = f.variables['time'] # A netCDF time variable object. + nctime = f.variables["time"] # A netCDF time variable object. idx = nc.date2index(date, nctime, select="exact") - idx = int(idx) + idx = int(idx) -# sameClone = True -# # check whether clone and input maps have the same attributes: -# if cloneMapFileName != None: -# # get the attributes of cloneMap -# #attributeClone = getMapAttributesALL(cloneMapFileName) -# #cellsizeClone = attributeClone['cellsize'] -# #rowsClone = attributeClone['rows'] -# #colsClone = attributeClone['cols'] -# #xULClone = attributeClone['xUL'] -# #yULClone = attributeClone['yUL'] -# attributeClone = getgridparams() -# cellsizeClone = attributeClone[2] -# rowsClone = attributeClone[4] -# colsClone = attributeClone[5] -# xULClone = attributeClone[0] - 0.5*cellsizeClone -# yULClone = attributeClone[1] + 0.5*cellsizeClone -# # get the attributes of input (netCDF) -# cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] -# cellsizeInput = float(cellsizeInput) -# rowsInput = len(f.variables['lat']) -# colsInput = len(f.variables['lon']) -# xULInput = f.variables['lon'][0]-0.5*cellsizeInput -# yULInput = f.variables['lat'][0]+0.5*cellsizeInput -# # check whether both maps have the same attributes -# if cellsizeClone != cellsizeInput: sameClone = False -# if rowsClone != rowsInput: sameClone = False -# if colsClone != colsInput: sameClone = False -# if xULClone != xULInput: sameClone = False -# if yULClone != yULInput: sameClone = False + # sameClone = True + # # check whether clone and input maps have the same attributes: + # if cloneMapFileName != None: + # # get the attributes of cloneMap + # #attributeClone = getMapAttributesALL(cloneMapFileName) + # #cellsizeClone = attributeClone['cellsize'] + # #rowsClone = attributeClone['rows'] + # #colsClone = attributeClone['cols'] + # #xULClone = attributeClone['xUL'] + # #yULClone = attributeClone['yUL'] + # attributeClone = getgridparams() + # cellsizeClone = attributeClone[2] + # rowsClone = attributeClone[4] + # colsClone = attributeClone[5] + # xULClone = attributeClone[0] - 0.5*cellsizeClone + # yULClone = attributeClone[1] + 0.5*cellsizeClone + # # get the attributes of input (netCDF) + # cellsizeInput = f.variables['lat'][0]- f.variables['lat'][1] + # cellsizeInput = float(cellsizeInput) + # rowsInput = len(f.variables['lat']) + # colsInput = len(f.variables['lon']) + # xULInput = f.variables['lon'][0]-0.5*cellsizeInput + # yULInput = f.variables['lat'][0]+0.5*cellsizeInput + # # check whether both maps have the same attributes + # if cellsizeClone != cellsizeInput: sameClone = False + # if rowsClone != rowsInput: sameClone = False + # if colsClone != colsInput: sameClone = False + # if xULClone != xULInput: sameClone = False + # if yULClone != yULInput: sameClone = False - cropData = f.variables[varName][int(idx),:,:] # still original data -# factor = 1 # needed in regridData2FinerGrid -# if sameClone == False: -# # crop to cloneMap: -# xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0]) -# xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) -# yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0]) -# yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) -# cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] -# factor = int(float(cellsizeInput)/float(cellsizeClone)) - + cropData = f.variables[varName][int(idx), :, :] # still original data + # factor = 1 # needed in regridData2FinerGrid + # if sameClone == False: + # # crop to cloneMap: + # xIdxSta = int(np.where(f.variables['lon'][:] == xULClone + 0.5*cellsizeInput)[0]) + # xIdxEnd = int(math.ceil(xIdxSta + colsClone /(cellsizeInput/cellsizeClone))) + # yIdxSta = int(np.where(f.variables['lat'][:] == yULClone - 0.5*cellsizeInput)[0]) + # yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone))) + # cropData = f.variables[varName][idx,yIdxSta:yIdxEnd,xIdxSta:xIdxEnd] + # factor = int(float(cellsizeInput)/float(cellsizeClone)) + # convert to PCR object and close f - outPCR = pcr.numpy2pcr(pcr.Scalar, cropData, \ -# regridData2FinerGrid(factor,cropData,MV), \ - float(f.variables[varName]._FillValue)) - f.close(); - f = None ; cropData = None + outPCR = pcr.numpy2pcr( + pcr.Scalar, + cropData, + # regridData2FinerGrid(factor,cropData,MV), \ + float(f.variables[varName]._FillValue), + ) + f.close() + f = None + cropData = None # PCRaster object - return (outPCR) - -def netcdf2PCRobj(ncFile,varName,dateInput): + return outPCR + + +def netcdf2PCRobj(ncFile, varName, dateInput): # EHS (04 APR 2013): To convert netCDF (tss) file to PCR file. # The cloneMap is globally defined (outside this method). - + # Get netCDF file and variable name: f = nc.Dataset(ncFile) varName = str(varName) # date date = dateInput - if isinstance(date, str) == True: date = \ - datetime.datetime.strptime(str(date),'%Y-%m-%d') - date = datetime.datetime(date.year,date.month,date.day) - + if isinstance(date, str) == True: + date = datetime.datetime.strptime(str(date), "%Y-%m-%d") + date = datetime.datetime(date.year, date.month, date.day) + # time index (in the netCDF file) - nctime = f.variables['time'] # A netCDF time variable object. - idx = nc.date2index(date, nctime, calendar=nctime.calendar, \ - select='exact') - + nctime = f.variables["time"] # A netCDF time variable object. + idx = nc.date2index(date, nctime, calendar=nctime.calendar, select="exact") + # convert to PCR object and close f - outPCR = pcr.numpy2pcr(pcr.Scalar,(f.variables[varName][idx].data), \ - float(f.variables[varName]._FillValue)) - f.close(); f = None ; del f + outPCR = pcr.numpy2pcr( + pcr.Scalar, + (f.variables[varName][idx].data), + float(f.variables[varName]._FillValue), + ) + f.close() + f = None + del f # PCRaster object - return (outPCR) + return outPCR + def makeDir(directoryName): try: os.makedirs(directoryName) except OSError: pass -def writePCRmapToDir(v,outFileName,outDir): + +def writePCRmapToDir(v, outFileName, outDir): # v: inputMapFileName or floating values # cloneMapFileName: If the inputMap and cloneMap have different clones, - # resampling will be done. Then, - fullFileName = getFullPath(outFileName,outDir) - logger.debug('Writing a pcraster map to : '+str(fullFileName)) - pcr.report(v,fullFileName) + # resampling will be done. Then, + fullFileName = getFullPath(outFileName, outDir) + logger.debug("Writing a pcraster map to : " + str(fullFileName)) + pcr.report(v, fullFileName) -def readPCRmapClone(v,cloneMapFileName,tmpDir,absolutePath=None,isLddMap=False,cover=None,isNomMap=False): - # v: inputMapFileName or floating values - # cloneMapFileName: If the inputMap and cloneMap have different clones, - # resampling will be done. - logger.debug('read file/values: '+str(v)) + +def readPCRmapClone( + v, + cloneMapFileName, + tmpDir, + absolutePath=None, + isLddMap=False, + cover=None, + isNomMap=False, +): + # v: inputMapFileName or floating values + # cloneMapFileName: If the inputMap and cloneMap have different clones, + # resampling will be done. + logger.debug("read file/values: " + str(v)) if v == "None": - #~ PCRmap = str("None") - PCRmap = None # 29 July: I made an experiment by changing the type of this object. - elif not re.match(r"[0-9.-]*$",v): - if absolutePath != None: v = getFullPath(v,absolutePath) + # ~ PCRmap = str("None") + PCRmap = ( + None + ) # 29 July: I made an experiment by changing the type of this object. + elif not re.match(r"[0-9.-]*$", v): + if absolutePath != None: + v = getFullPath(v, absolutePath) # print(v) - sameClone = True #isSameClone(v,cloneMapFileName) #for wflow CloneMap and inputMap should be the same + sameClone = ( + True + ) # isSameClone(v,cloneMapFileName) #for wflow CloneMap and inputMap should be the same if sameClone == True: PCRmap = pcr.readmap(v) else: # resample using GDAL: - output = tmpDir+'temp.map' - warp = gdalwarpPCR(v,output,cloneMapFileName,tmpDir,isLddMap,isNomMap) + output = tmpDir + "temp.map" + warp = gdalwarpPCR(v, output, cloneMapFileName, tmpDir, isLddMap, isNomMap) # read from temporary file and delete the temporary file: PCRmap = pcr.readmap(output) - if isLddMap == True: PCRmap = pcr.ifthen(pcr.scalar(PCRmap) < 10., PCRmap) - if isLddMap == True: PCRmap = pcr.ldd(PCRmap) - if isNomMap == True: PCRmap = pcr.ifthen(pcr.scalar(PCRmap) > 0., PCRmap) - if isNomMap == True: PCRmap = pcr.nominal(PCRmap) + if isLddMap == True: + PCRmap = pcr.ifthen(pcr.scalar(PCRmap) < 10., PCRmap) + if isLddMap == True: + PCRmap = pcr.ldd(PCRmap) + if isNomMap == True: + PCRmap = pcr.ifthen(pcr.scalar(PCRmap) > 0., PCRmap) + if isNomMap == True: + PCRmap = pcr.nominal(PCRmap) if os.path.isdir(tmpDir): shutil.rmtree(tmpDir) os.makedirs(tmpDir) else: PCRmap = pcr.spatial(pcr.scalar(float(v))) if cover != None: PCRmap = pcr.cover(PCRmap, cover) - co = None; cOut = None; err = None; warp = None - del co; del cOut; del err; del warp - stdout = None; del stdout - stderr = None; del stderr - return PCRmap + co = None + cOut = None + err = None + warp = None + del co + del cOut + del err + del warp + stdout = None + del stdout + stderr = None + del stderr + return PCRmap + def readPCRmap(v): - # v : fileName or floating values + # v : fileName or floating values if not re.match(r"[0-9.-]*$", v): PCRmap = pcr.readmap(v) else: PCRmap = pcr.scalar(float(v)) - return PCRmap + return PCRmap -def isSameClone(inputMapFileName,cloneMapFileName): + +def isSameClone(inputMapFileName, cloneMapFileName): # reading inputMap: attributeInput = getMapAttributesALL(inputMapFileName) - cellsizeInput = attributeInput['cellsize'] - rowsInput = attributeInput['rows'] - colsInput = attributeInput['cols'] - xULInput = attributeInput['xUL'] - yULInput = attributeInput['yUL'] + cellsizeInput = attributeInput["cellsize"] + rowsInput = attributeInput["rows"] + colsInput = attributeInput["cols"] + xULInput = attributeInput["xUL"] + yULInput = attributeInput["yUL"] # reading cloneMap: attributeClone = getMapAttributesALL(cloneMapFileName) - cellsizeClone = attributeClone['cellsize'] - rowsClone = attributeClone['rows'] - colsClone = attributeClone['cols'] - xULClone = attributeClone['xUL'] - yULClone = attributeClone['yUL'] - # check whether both maps have the same attributes? + cellsizeClone = attributeClone["cellsize"] + rowsClone = attributeClone["rows"] + colsClone = attributeClone["cols"] + xULClone = attributeClone["xUL"] + yULClone = attributeClone["yUL"] + # check whether both maps have the same attributes? sameClone = True - if cellsizeClone != cellsizeInput: sameClone = False - if rowsClone != rowsInput: sameClone = False - if colsClone != colsInput: sameClone = False - if xULClone != xULInput: sameClone = False - if yULClone != yULInput: sameClone = False + if cellsizeClone != cellsizeInput: + sameClone = False + if rowsClone != rowsInput: + sameClone = False + if colsClone != colsInput: + sameClone = False + if xULClone != xULInput: + sameClone = False + if yULClone != yULInput: + sameClone = False return sameClone -def gdalwarpPCR(input,output,cloneOut,tmpDir,isLddMap=False,isNominalMap=False): + +def gdalwarpPCR(input, output, cloneOut, tmpDir, isLddMap=False, isNominalMap=False): # 19 Mar 2013 created by Edwin H. Sutanudjaja # all input maps must be in PCRaster maps - # + # # remove temporary files: - co = 'rm '+str(tmpDir)+'*.*' - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - # + co = "rm " + str(tmpDir) + "*.*" + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open(os.devnull), shell=True + ).communicate() + # # converting files to tif: - co = 'gdal_translate -ot Float64 '+str(input)+' '+str(tmpDir)+'tmp_inp.tif' - if isLddMap == True: co = 'gdal_translate -ot Int32 '+str(input)+' '+str(tmpDir)+'tmp_inp.tif' - if isNominalMap == True: co = 'gdal_translate -ot Int32 '+str(input)+' '+str(tmpDir)+'tmp_inp.tif' - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - # + co = "gdal_translate -ot Float64 " + str(input) + " " + str(tmpDir) + "tmp_inp.tif" + if isLddMap == True: + co = ( + "gdal_translate -ot Int32 " + str(input) + " " + str(tmpDir) + "tmp_inp.tif" + ) + if isNominalMap == True: + co = ( + "gdal_translate -ot Int32 " + str(input) + " " + str(tmpDir) + "tmp_inp.tif" + ) + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open(os.devnull), shell=True + ).communicate() + # # get the attributes of PCRaster map: cloneAtt = getMapAttributesALL(cloneOut) - xmin = cloneAtt['xUL'] - ymin = cloneAtt['yUL'] - cloneAtt['rows']*cloneAtt['cellsize'] - xmax = cloneAtt['xUL'] + cloneAtt['cols']*cloneAtt['cellsize'] - ymax = cloneAtt['yUL'] - xres = cloneAtt['cellsize'] - yres = cloneAtt['cellsize'] - te = '-te '+str(xmin)+' '+str(ymin)+' '+str(xmax)+' '+str(ymax)+' ' - tr = '-tr '+str(xres)+' '+str(yres)+' ' - co = 'gdalwarp '+te+tr+ \ - ' -srcnodata -3.4028234663852886e+38 -dstnodata mv '+ \ - str(tmpDir)+'tmp_inp.tif '+ \ - str(tmpDir)+'tmp_out.tif' - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - # - co = 'gdal_translate -of PCRaster '+ \ - str(tmpDir)+'tmp_out.tif '+str(output) - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - # - co = 'mapattr -c '+str(cloneOut)+' '+str(output) - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - # - #~ co = 'aguila '+str(output) - #~ print(co) - #~ cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - # - co = 'rm '+str(tmpDir)+'tmp*.*' - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() - co = None; cOut = None; err = None - del co; del cOut; del err - stdout = None; del stdout - stderr = None; del stderr - n = gc.collect() ; del gc.garbage[:] ; n = None ; del n + xmin = cloneAtt["xUL"] + ymin = cloneAtt["yUL"] - cloneAtt["rows"] * cloneAtt["cellsize"] + xmax = cloneAtt["xUL"] + cloneAtt["cols"] * cloneAtt["cellsize"] + ymax = cloneAtt["yUL"] + xres = cloneAtt["cellsize"] + yres = cloneAtt["cellsize"] + te = "-te " + str(xmin) + " " + str(ymin) + " " + str(xmax) + " " + str(ymax) + " " + tr = "-tr " + str(xres) + " " + str(yres) + " " + co = ( + "gdalwarp " + + te + + tr + + " -srcnodata -3.4028234663852886e+38 -dstnodata mv " + + str(tmpDir) + + "tmp_inp.tif " + + str(tmpDir) + + "tmp_out.tif" + ) + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open(os.devnull), shell=True + ).communicate() + # + co = "gdal_translate -of PCRaster " + str(tmpDir) + "tmp_out.tif " + str(output) + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open(os.devnull), shell=True + ).communicate() + # + co = "mapattr -c " + str(cloneOut) + " " + str(output) + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open(os.devnull), shell=True + ).communicate() + # + # ~ co = 'aguila '+str(output) + # ~ print(co) + # ~ cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open(os.devnull),shell=True).communicate() + # + co = "rm " + str(tmpDir) + "tmp*.*" + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open(os.devnull), shell=True + ).communicate() + co = None + cOut = None + err = None + del co + del cOut + del err + stdout = None + del stdout + stderr = None + del stderr + n = gc.collect() + del gc.garbage[:] + n = None + del n -def getFullPath(inputPath,absolutePath,completeFileName = True): + +def getFullPath(inputPath, absolutePath, completeFileName=True): # 19 Mar 2013 created by Edwin H. Sutanudjaja # Function: to get the full absolute path of a folder or a file - + # replace all \ with / inputPath = str(inputPath).replace("\\", "/") absolutePath = str(absolutePath).replace("\\", "/") - + # tuple of suffixes (extensions) that can be used: - suffix = ('/','_','.nc4','.map','.nc','.dat','.txt','.asc','.ldd','.tbl',\ - '.001','.002','.003','.004','.005','.006',\ - '.007','.008','.009','.010','.011','.012') - - if inputPath.startswith('/') or str(inputPath)[1] == ":": + suffix = ( + "/", + "_", + ".nc4", + ".map", + ".nc", + ".dat", + ".txt", + ".asc", + ".ldd", + ".tbl", + ".001", + ".002", + ".003", + ".004", + ".005", + ".006", + ".007", + ".008", + ".009", + ".010", + ".011", + ".012", + ) + + if inputPath.startswith("/") or str(inputPath)[1] == ":": fullPath = str(inputPath) else: - if absolutePath.endswith('/'): + if absolutePath.endswith("/"): absolutePath = str(absolutePath) else: - absolutePath = str(absolutePath)+'/' - fullPath = str(absolutePath)+str(inputPath) - + absolutePath = str(absolutePath) + "/" + fullPath = str(absolutePath) + str(inputPath) + if completeFileName: - if fullPath.endswith(suffix): + if fullPath.endswith(suffix): fullPath = str(fullPath) - else: - fullPath = str(fullPath)+'/' + else: + fullPath = str(fullPath) + "/" - return fullPath + return fullPath -def findISIFileName(year,model,rcp,prefix,var): - histYears = [1951,1961,1971,1981,1991,2001] - sYears = [2011,2021,2031,2041,2051,2061,2071,2081,2091] + +def findISIFileName(year, model, rcp, prefix, var): + histYears = [1951, 1961, 1971, 1981, 1991, 2001] + sYears = [2011, 2021, 2031, 2041, 2051, 2061, 2071, 2081, 2091] rcpStr = rcp if year >= sYears[0]: sYear = [i for i in range(len(sYears)) if year >= sYears[i]] - sY = sYears[sYear[-1]] - + sY = sYears[sYear[-1]] + elif year < histYears[-1]: - - sYear = [i for i in range(len(histYears)) if year >= histYears[i] ] - sY = histYears[sYear[-1]] - + + sYear = [i for i in range(len(histYears)) if year >= histYears[i]] + sY = histYears[sYear[-1]] + if year >= histYears[-1] and year < sYears[0]: - - if model == 'HadGEM2-ES': + + if model == "HadGEM2-ES": if year < 2005: - rcpStr = 'historical' + rcpStr = "historical" sY = 2001 eY = 2004 else: rcpStr = rcp sY = 2005 eY = 2010 - if model == 'IPSL-CM5A-LR' or model == 'GFDL-ESM2M': + if model == "IPSL-CM5A-LR" or model == "GFDL-ESM2M": if year < 2006: - rcpStr = 'historical' + rcpStr = "historical" sY = 2001 eY = 2005 else: rcpStr = rcp sY = 2006 eY = 2010 - - else: + + else: eY = sY + 9 if sY == 2091: - eY = 2099 - if model == 'HadGEM2-ES': + eY = 2099 + if model == "HadGEM2-ES": if year < 2005: - rcpStr = 'historical' - if model == 'IPSL-CM5A-LR' or model == 'GFDL-ESM2M': + rcpStr = "historical" + if model == "IPSL-CM5A-LR" or model == "GFDL-ESM2M": if year < 2006: - rcpStr = 'historical' - #print year,sY,eY - return "%s_%s_%s_%s_%i-%i.nc" %(var,prefix,model.lower(),rcpStr,sY,eY) + rcpStr = "historical" + # print year,sY,eY + return "%s_%s_%s_%s_%i-%i.nc" % (var, prefix, model.lower(), rcpStr, sY, eY) - + def get_random_word(wordLen): - word = '' + word = "" for i in range(wordLen): - word += random.choice('ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789') + word += random.choice( + "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789" + ) return word - + + def isLastDayOfMonth(date): - if (date + datetime.timedelta(days=1 )).day == 1: + if (date + datetime.timedelta(days=1)).day == 1: return True else: return False -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 ? ") +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; err = None - del co; del err - n = gc.collect() ; del gc.garbage[:] ; n = None ; del n - if attribute == 'cellsize': + # print cOut.split() + co = None + err = None + del co + del err + n = gc.collect() + del gc.garbage[:] + n = None + del n + if attribute == "cellsize": cellsize = float(cOut.split()[7]) - if arcDegree == True: cellsize = round(cellsize * 360000.)/360000. - return cellsize - if attribute == 'rows': + if arcDegree == True: + cellsize = round(cellsize * 360000.) / 360000. + return cellsize + if attribute == "rows": return int(cOut.split()[3]) - #return float(cOut.split()[3]) - if attribute == 'cols': + # return float(cOut.split()[3]) + if attribute == "cols": return int(cOut.split()[5]) - #return float(cOut.split()[5]) - if attribute == 'xUL': + # return float(cOut.split()[5]) + if attribute == "xUL": return float(cOut.split()[17]) - if attribute == 'yUL': + if attribute == "yUL": return float(cOut.split()[19]) - + + def getMapTotal(mapFile): - ''' outputs the sum of all values in a map file ''' + """ outputs the sum of all values in a map file """ - total, valid = pcr.cellvalue(pcr.maptotal(mapFile),1) + total, valid = pcr.cellvalue(pcr.maptotal(mapFile), 1) return total + def getMapTotalHighPrecisionButOnlyForPositiveValues_NEEDMORETEST(mapFile): - ''' outputs the sum of all values in a map file ''' + """ outputs the sum of all values in a map file """ # STILL UNDER DEVELOPMENT - NOT FULLY TESTED - + # input map - note that all values must be positive remainingMapValue = pcr.max(0.0, mapFile) - + # loop from biggest values - min_power_number = 0 # The minimum value is zero. + min_power_number = 0 # The minimum value is zero. max_power_number = int(pcr.mapmaximum(pcr.log10(remainingMapValue))) + 1 step = 1 total_map_for_every_power_number = {} for power_number in range(max_power_number, min_power_number - step, -step): - - # cell value in this loop - currentCellValue = pcr.rounddown(remainingMapValue * pcr.scalar(10.**(power_number))) / pcr.scalar(10.**(power_number)) - if power_number == min_power_number: currentCellValue = remainingMapValue - + + # cell value in this loop + currentCellValue = pcr.rounddown( + remainingMapValue * pcr.scalar(10. ** (power_number)) + ) / pcr.scalar(10. ** (power_number)) + if power_number == min_power_number: + currentCellValue = remainingMapValue + # map total in this loop total_in_this_loop, valid = pcr.cellvalue(pcr.maptotal(currentCellValue), 1) total_map_for_every_power_number[str(power_number)] = total_in_this_loop - - # remaining map value + + # remaining map value remainingMapValue = pcr.max(0.0, remainingMapValue - currentCellValue) - + # sum from the smallest values (minimizing numerical errors) total = pcr.scalar(0.0) for power_number in range(min_power_number, max_power_number + step, step): total += total_map_for_every_power_number[str(power_number)] return total + def get_rowColAboveThreshold(map, threshold): npMap = pcr.pcr2numpy(map, -9999) (nr, nc) = np.shape(npMap) @@ -1154,193 +1546,240 @@ if npMap[r, c] != -9999: if np.abs(npMap[r, c]) > threshold: - return (r, c) def getLastDayOfMonth(date): - ''' returns the last day of the month for a given date ''' + """ returns the last day of the month for a given date """ if date.month == 12: return date.replace(day=31) return date.replace(month=date.month + 1, day=1) - datetime.timedelta(days=1) - -def getMinMaxMean(mapFile,ignoreEmptyMap=False): - mn = pcr.cellvalue(pcr.mapminimum(mapFile),1)[0] - mx = pcr.cellvalue(pcr.mapmaximum(mapFile),1)[0] - nrValues = pcr.cellvalue(pcr.maptotal(pcr.scalar(pcr.defined(mapFile))), 1 )[0] #/ getNumNonMissingValues(mapFile) - if nrValues == 0.0 and ignoreEmptyMap: - return 0.0,0.0,0.0 +def getMinMaxMean(mapFile, ignoreEmptyMap=False): + mn = pcr.cellvalue(pcr.mapminimum(mapFile), 1)[0] + mx = pcr.cellvalue(pcr.mapmaximum(mapFile), 1)[0] + nrValues = pcr.cellvalue(pcr.maptotal(pcr.scalar(pcr.defined(mapFile))), 1)[ + 0 + ] # / getNumNonMissingValues(mapFile) + if nrValues == 0.0 and ignoreEmptyMap: + return 0.0, 0.0, 0.0 else: - return mn,mx,(getMapTotal(mapFile) / nrValues) + return mn, mx, (getMapTotal(mapFile) / nrValues) + def getMapVolume(mapFile, cellareaFile): - ''' returns the sum of all grid cell values ''' + """ returns the sum of all grid cell values """ volume = mapFile * cellareaFile - return (getMapTotal(volume) / 1) + return getMapTotal(volume) / 1 + 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) + +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 range(nrRows): - row,col= coordinates[iCnt,:] - if row != MV and col != MV: - x[iCnt]= tmpIDArray[row,col] + row, col = coordinates[iCnt, :] + if row != MV and col != MV: + x[iCnt] = tmpIDArray[row, col] return x -def returnMapValue(pcrX,x,coord): - #-retrieves value from an array and update values in the map + +def returnMapValue(pcrX, x, coord): + # -retrieves value from an array and update values in the map if x.ndim == 1: - nrRows= 1 + nrRows = 1 - tempIDArray= pcr.pcr2numpy(pcrX,MV) - #print tempIDArray - temporary= tempIDArray - nrRows= coord.shape[0] + 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) + 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 - + + def getQAtBasinMouths(discharge, basinMouth): - temp = pcr.ifthenelse(basinMouth != 0 , discharge * secondsPerDay(),0.) - pcr.report(temp,"temp.map") - return (getMapTotal(temp) / 1e9) + temp = pcr.ifthenelse(basinMouth != 0, discharge * secondsPerDay(), 0.) + pcr.report(temp, "temp.map") + return getMapTotal(temp) / 1e9 -def regridMapFile2FinerGrid (rescaleFac,coarse): - if rescaleFac ==1: + +def regridMapFile2FinerGrid(rescaleFac, coarse): + if rescaleFac == 1: return coarse - return pcr.numpy2pcr(pcr.Scalar, regridData2FinerGrid(rescaleFac,pcr.pcr2numpy(coarse,MV),MV),MV) - -def regridData2FinerGrid(rescaleFac,coarse,MV): - if rescaleFac ==1: + return pcr.numpy2pcr( + pcr.Scalar, regridData2FinerGrid(rescaleFac, pcr.pcr2numpy(coarse, MV), MV), MV + ) + + +def regridData2FinerGrid(rescaleFac, coarse, MV): + if rescaleFac == 1: return coarse - nr,nc = np.shape(coarse) - - fine= np.zeros(nr*nc*rescaleFac*rescaleFac).reshape(nr*rescaleFac,nc*rescaleFac) + MV - - + nr, nc = np.shape(coarse) + + fine = ( + np.zeros(nr * nc * rescaleFac * rescaleFac).reshape( + nr * rescaleFac, nc * rescaleFac + ) + + MV + ) + ii = -1 - nrF,ncF = np.shape(fine) - for i in range(0 , nrF): - if i % rescaleFac == 0: - ii += 1 - fine [i,:] = coarse[ii,:].repeat(rescaleFac) + nrF, ncF = np.shape(fine) + for i in range(0, nrF): + if i % rescaleFac == 0: + ii += 1 + fine[i, :] = coarse[ii, :].repeat(rescaleFac) - nr = None; nc = None - del nr; del nc - nrF = None; ncF = None - del nrF; del ncF - n = gc.collect() ; del gc.garbage[:] ; n = None ; del n + nr = None + nc = None + del nr + del nc + nrF = None + ncF = None + del nrF + del ncF + n = gc.collect() + del gc.garbage[:] + n = None + del n return fine -def regridToCoarse(fine,fac,mode,missValue): - nr,nc = np.shape(fine) - coarse = np.zeros(nr/fac * nc / fac).reshape(nr/fac,nc/fac) + MV - nr,nc = np.shape(coarse) - for r in range(0,nr): - for c in range(0,nc): - ar = fine[r * fac : fac * (r+1),c * fac: fac * (c+1)] - m = np.ma.masked_values(ar,missValue) + +def regridToCoarse(fine, fac, mode, missValue): + nr, nc = np.shape(fine) + coarse = np.zeros(nr / fac * nc / fac).reshape(nr / fac, nc / fac) + MV + nr, nc = np.shape(coarse) + for r in range(0, nr): + for c in range(0, nc): + ar = fine[r * fac : fac * (r + 1), c * fac : fac * (c + 1)] + m = np.ma.masked_values(ar, missValue) if ma.count(m) == 0: - coarse[r,c] = MV + coarse[r, c] = MV else: - if mode == 'average': - coarse [r,c] = ma.average(m) - elif mode == 'median': - coarse [r,c] = ma.median(m) - elif mode == 'sum': - coarse [r,c] = ma.sum(m) - elif mode =='min': - coarse [r,c] = ma.min(m) - elif mode == 'max': - coarse [r,c] = ma.max(m) - return coarse - - -def waterBalanceCheck(fluxesIn,fluxesOut,preStorages,endStorages,processName,PrintOnlyErrors,dateStr,threshold=1e-5,landmask=None): + if mode == "average": + coarse[r, c] = ma.average(m) + elif mode == "median": + coarse[r, c] = ma.median(m) + elif mode == "sum": + coarse[r, c] = ma.sum(m) + elif mode == "min": + coarse[r, c] = ma.min(m) + elif mode == "max": + coarse[r, c] = ma.max(m) + return coarse + + +def waterBalanceCheck( + fluxesIn, + fluxesOut, + preStorages, + endStorages, + processName, + PrintOnlyErrors, + dateStr, + threshold=1e-5, + landmask=None, +): """ Returns the water balance for a list of input, output, and storage map files """ # modified by Edwin (22 Apr 2013) - inMap = pcr.spatial(pcr.scalar(0.0)) - outMap = pcr.spatial(pcr.scalar(0.0)) - dsMap = pcr.spatial(pcr.scalar(0.0)) - + inMap = pcr.spatial(pcr.scalar(0.0)) + outMap = pcr.spatial(pcr.scalar(0.0)) + dsMap = pcr.spatial(pcr.scalar(0.0)) + for fluxIn in fluxesIn: - inMap += fluxIn + inMap += fluxIn for fluxOut in fluxesOut: - outMap += fluxOut + outMap += fluxOut for preStorage in preStorages: - dsMap += preStorage + dsMap += preStorage for endStorage in endStorages: - dsMap -= endStorage + dsMap -= endStorage - a,b,c = getMinMaxMean(inMap + dsMap- outMap) + a, b, c = getMinMaxMean(inMap + dsMap - outMap) if abs(a) > threshold or abs(b) > threshold: - if PrintOnlyErrors: - - msg = "\n" + if PrintOnlyErrors: + + msg = "\n" msg += "\n" - msg = "\n" + msg = "\n" msg += "\n" - msg += "##############################################################################################################################################\n" - msg += "WARNING !!!!!!!! Water Balance Error %s Min %f Max %f Mean %f" %(processName,a,b,c) + msg += ( + "##############################################################################################################################################\n" + ) + msg += "WARNING !!!!!!!! Water Balance Error %s Min %f Max %f Mean %f" % ( + processName, + a, + b, + c, + ) msg += "\n" - msg += "##############################################################################################################################################\n" + msg += ( + "##############################################################################################################################################\n" + ) msg += "\n" msg += "\n" msg += "\n" - + logger.error(msg) - #~ pcr.report(inMap + dsMap - outMap,"wb.map") - #~ os.system("aguila wb.map") - - #~ # for debugging: - #~ error = inMap + dsMap- outMap - #~ os.system('rm error.map') - #~ pcr.report(error,"error.map") - #~ os.system('aguila error.map') - #~ os.system('rm error.map') - - #~ wb = inMap + dsMap - outMap - #~ maxWBError = pcr.cellvalue(pcr.mapmaximum(pcr.abs(wb)), 1, 1)[0] - #~ #return wb + # ~ pcr.report(inMap + dsMap - outMap,"wb.map") + # ~ os.system("aguila wb.map") + # ~ # for debugging: + # ~ error = inMap + dsMap- outMap + # ~ os.system('rm error.map') + # ~ pcr.report(error,"error.map") + # ~ os.system('aguila error.map') + # ~ os.system('rm error.map') + # ~ wb = inMap + dsMap - outMap + # ~ maxWBError = pcr.cellvalue(pcr.mapmaximum(pcr.abs(wb)), 1, 1)[0] + # ~ #return wb -def waterBalance( fluxesIn, fluxesOut, deltaStorages, processName, PrintOnlyErrors, dateStr,threshold=1e-5): +def waterBalance( + fluxesIn, + fluxesOut, + deltaStorages, + processName, + PrintOnlyErrors, + dateStr, + threshold=1e-5, +): """ Returns the water balance for a list of input, output, and storage map files and """ inMap = pcr.spatial(pcr.scalar(0.0)) @@ -1359,249 +1798,315 @@ deltaS += getMapTotal(deltaStorage) dsMap += deltaStorage - #if PrintOnlyErrors: - a,b,c = getMinMaxMean(inMap + dsMap- outMap) + # if PrintOnlyErrors: + a, b, c = getMinMaxMean(inMap + dsMap - outMap) # 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) - #else: + # else: # print "Water balance for %s: on %s in = %f\tout=%f\tdeltaS=%f\tBalance=%f" \ # %(processName,dateStr,inflow,outflow,deltaS,inflow + deltaS - outflow) wb = inMap + dsMap - outMap maxWBError = pcr.cellvalue(pcr.mapmaximum(pcr.abs(wb)), 1, 1)[0] - #if maxWBError > 0.001 / 1000: - #row = 0 - #col = 0 - #cellID = 1 - #troubleCell = 0 + # if maxWBError > 0.001 / 1000: + # row = 0 + # col = 0 + # cellID = 1 + # troubleCell = 0 - #print "Water balance for %s on %s: %f mm !!! " %(processName,dateStr,maxWBError * 1000) - #pcr.report(wb,"%s-WaterBalanceError-%s" %(processName,dateStr)) + # print "Water balance for %s on %s: %f mm !!! " %(processName,dateStr,maxWBError * 1000) + # pcr.report(wb,"%s-WaterBalanceError-%s" %(processName,dateStr)) - #npWBMError = pcr2numpy(wb, -9999) - #(nr, nc) = np.shape(npWBMError) - #for r in range(0, nr): - #for c in range(0, nc): + # npWBMError = pcr2numpy(wb, -9999) + # (nr, nc) = np.shape(npWBMError) + # for r in range(0, nr): + # for c in range(0, nc): - ## print r,c + ## print r,c - #if npWBMError[r, c] != -9999.0: - #val = npWBMError[r, c] - #if math.fabs(val) > 0.0001 / 1000: + # if npWBMError[r, c] != -9999.0: + # val = npWBMError[r, c] + # if math.fabs(val) > 0.0001 / 1000: - ## print npWBMError[r,c] + ## print npWBMError[r,c] - #row = r - #col = c - #troubleCell = cellID - #cellID += 1 - #print 'Water balance for %s on %s: %f mm row %i col %i cellID %i!!! ' % ( - #processName, - #dateStr, - #maxWBError * 1000, - #row, - #col, - #troubleCell, - #) + # row = r + # col = c + # troubleCell = cellID + # cellID += 1 + # print 'Water balance for %s on %s: %f mm row %i col %i cellID %i!!! ' % ( + # processName, + # dateStr, + # maxWBError * 1000, + # row, + # col, + # troubleCell, + # ) return inMap + dsMap - outMap +def waterAbstractionAndAllocationHighPrecision_NEEDMORETEST( + water_demand_volume, + available_water_volume, + allocation_zones, + zone_area=None, + debug_water_balance=True, + extra_info_for_water_balance_reporting="", +): -def waterAbstractionAndAllocationHighPrecision_NEEDMORETEST(water_demand_volume, \ - available_water_volume, \ - allocation_zones,\ - zone_area = None, - debug_water_balance = True,\ - extra_info_for_water_balance_reporting = ""): - # STILL UNDER DEVELOPMENT - NOT FULLY TESTED - + logger.debug("Allocation of abstraction. - using high precision option") - + # demand volume in each cell (unit: m3) remainingcellVolDemand = pcr.max(0.0, water_demand_volume) - + # available water volume in each cell - remainingCellAvlWater = pcr.max(0.0, available_water_volume) + remainingCellAvlWater = pcr.max(0.0, available_water_volume) # loop from biggest values of cellAvlWater - min_power_number = 0 # The minimum value is zero. + min_power_number = 0 # The minimum value is zero. max_power_number = int(pcr.mapmaximum(pcr.log10(remainingCellAvlWater))) + 1 step = 1 cell_abstrac_for_every_power_number = {} cell_allocat_for_every_power_number = {} for power_number in range(max_power_number, min_power_number - step, -step): - - logger.debug("Allocation of abstraction. - using high precision option - loop power number: " + str(power_number)) + logger.debug( + "Allocation of abstraction. - using high precision option - loop power number: " + + str(power_number) + ) - # cell available water in this loop - cellAvlWater = pcr.rounddown(remainingCellAvlWater * pcr.scalar(10.**(power_number))) / pcr.scalar(10.**(power_number)) - if power_number == min_power_number: cellAvlWater = pcr.max(0.0, remainingCellAvlWater) - + # cell available water in this loop + cellAvlWater = pcr.rounddown( + remainingCellAvlWater * pcr.scalar(10. ** (power_number)) + ) / pcr.scalar(10. ** (power_number)) + if power_number == min_power_number: + cellAvlWater = pcr.max(0.0, remainingCellAvlWater) + # zonal available water in this loop zoneAvlWater = pcr.areatotal(cellAvlWater, allocation_zones) # total actual water abstraction volume in each zone/segment (unit: m3) # - limited to available water - zoneVolDemand = pcr.areatotal(remainingcellVolDemand, allocation_zones) + zoneVolDemand = pcr.areatotal(remainingcellVolDemand, allocation_zones) zoneAbstraction = pcr.min(zoneAvlWater, zoneVolDemand) - + # actual water abstraction volume in each cell (unit: m3) - cellAbstraction = getValDivZero(\ - cellAvlWater, zoneAvlWater, smallNumber) * zoneAbstraction - cellAbstraction = pcr.min(cellAbstraction, cellAvlWater) - + cellAbstraction = ( + getValDivZero(cellAvlWater, zoneAvlWater, smallNumber) * zoneAbstraction + ) + cellAbstraction = pcr.min(cellAbstraction, cellAvlWater) + # allocation water to meet water demand (unit: m3) - cellAllocation = getValDivZero(\ - remainingcellVolDemand, zoneVolDemand, smallNumber) * zoneAbstraction - + cellAllocation = ( + getValDivZero(remainingcellVolDemand, zoneVolDemand, smallNumber) + * zoneAbstraction + ) + # 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, 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, + ) # actual water abstraction and allocation in this current loop (power number) cell_abstrac_for_every_power_number[str(power_number)] = cellAbstraction cell_allocat_for_every_power_number[str(power_number)] = cellAllocation - - # remaining cell available water and demand - remainingCellAvlWater = pcr.max(0.0, remainingCellAvlWater - cellAbstraction) - remainingcellVolDemand = pcr.max(0.0, remainingcellVolDemand - cellAllocation ) - + + # remaining cell available water and demand + remainingCellAvlWater = pcr.max(0.0, remainingCellAvlWater - cellAbstraction) + remainingcellVolDemand = pcr.max(0.0, remainingcellVolDemand - cellAllocation) + # sum from the smallest values (minimizing numerical errors) sumCellAbstraction = pcr.scalar(0.0) - sumCellAllocation = pcr.scalar(0.0) + sumCellAllocation = pcr.scalar(0.0) for power_number in range(min_power_number, max_power_number + step, step): sumCellAbstraction += cell_abstrac_for_every_power_number[str(power_number)] - sumCellAllocation += cell_allocat_for_every_power_number[str(power_number)] - + 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, 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, + ) + return sumCellAbstraction, sumCellAllocation -def waterAbstractionAndAllocation(water_demand_volume,available_water_volume,allocation_zones,\ - zone_area = None, - high_volume_treshold = 1000000., - debug_water_balance = True,\ - extra_info_for_water_balance_reporting = "", - landmask = None, - ignore_small_values = False): +def waterAbstractionAndAllocation( + water_demand_volume, + available_water_volume, + allocation_zones, + zone_area=None, + high_volume_treshold=1000000., + debug_water_balance=True, + extra_info_for_water_balance_reporting="", + landmask=None, + ignore_small_values=False, +): + logger.debug("Allocation of abstraction.") - + # demand volume in each cell (unit: m3) cellVolDemand = pcr.max(0.0, water_demand_volume) if not isinstance(landmask, type(None)): cellVolDemand = pcr.ifthen(landmask, pcr.cover(cellVolDemand, 0.0)) - if ignore_small_values: # ignore small values to avoid runding error + if ignore_small_values: # ignore small values to avoid runding error cellVolDemand = pcr.rounddown(pcr.max(0.0, water_demand_volume)) else: cellVolDemand = pcr.max(0.0, water_demand_volume) - + # total demand volume in each zone/segment (unit: m3) zoneVolDemand = pcr.areatotal(cellVolDemand, allocation_zones) - + # total available water volume in each cell cellAvlWater = pcr.max(0.0, available_water_volume) if not isinstance(landmask, type(None)): cellAvlWater = pcr.ifthen(landmask, pcr.cover(cellAvlWater, 0.0)) - if ignore_small_values: # ignore small values to avoid runding error + 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, type(None)): # 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)), pcr.boolean(0)) - zoneAvlWater = pcr.areatotal( - pcr.ifthenelse(mask, 0.0, cellAvlWater), allocation_zones) - zoneAvlWater += pcr.areatotal( - pcr.ifthenelse(mask, cellAvlWater, 0.0), allocation_zones) + mask = pcr.cover( + pcr.ifthen(cellAvlWater > high_volume_treshold, pcr.boolean(1)), + pcr.boolean(0), + ) + zoneAvlWater = pcr.areatotal( + pcr.ifthenelse(mask, 0.0, cellAvlWater), allocation_zones + ) + zoneAvlWater += pcr.areatotal( + pcr.ifthenelse(mask, cellAvlWater, 0.0), allocation_zones + ) else: - zoneAvlWater = pcr.areatotal(cellAvlWater, allocation_zones) - + zoneAvlWater = pcr.areatotal(cellAvlWater, allocation_zones) + # total actual water abstraction volume in each zone/segment (unit: m3) # - limited to available water zoneAbstraction = pcr.min(zoneAvlWater, zoneVolDemand) - + # actual water abstraction volume in each cell (unit: m3) - cellAbstraction = getValDivZero(\ - cellAvlWater, zoneAvlWater, smallNumber)*zoneAbstraction - cellAbstraction = pcr.min(cellAbstraction, cellAvlWater) - if ignore_small_values: # ignore small values to avoid runding error + cellAbstraction = ( + getValDivZero(cellAvlWater, zoneAvlWater, smallNumber) * zoneAbstraction + ) + 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, type(None)): # 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)), pcr.boolean(0)) - zoneAbstraction = pcr.areatotal( - pcr.ifthenelse(mask, 0.0, cellAbstraction), allocation_zones) - zoneAbstraction += pcr.areatotal( - pcr.ifthenelse(mask, cellAbstraction, 0.0), allocation_zones) + mask = pcr.cover( + pcr.ifthen(cellAbstraction > high_volume_treshold, pcr.boolean(1)), + pcr.boolean(0), + ) + zoneAbstraction = pcr.areatotal( + pcr.ifthenelse(mask, 0.0, cellAbstraction), allocation_zones + ) + zoneAbstraction += pcr.areatotal( + pcr.ifthenelse(mask, cellAbstraction, 0.0), allocation_zones + ) else: - zoneAbstraction = pcr.areatotal(cellAbstraction, allocation_zones) - + zoneAbstraction = pcr.areatotal(cellAbstraction, allocation_zones) + # allocation water to meet water demand (unit: m3) - cellAllocation = getValDivZero(\ - cellVolDemand, zoneVolDemand, smallNumber)*zoneAbstraction - + cellAllocation = ( + getValDivZero(cellVolDemand, zoneVolDemand, smallNumber) * zoneAbstraction + ) + # extraAbstraction to minimize numerical errors: - zoneDeficitAbstraction = pcr.max(0.0,\ - pcr.areatotal(cellAllocation , allocation_zones) -\ - pcr.areatotal(cellAbstraction, allocation_zones)) + zoneDeficitAbstraction = pcr.max( + 0.0, + pcr.areatotal(cellAllocation, allocation_zones) + - pcr.areatotal(cellAbstraction, allocation_zones), + ) remainingCellAvlWater = pcr.max(0.0, cellAvlWater - cellAbstraction) - cellAbstraction += zoneDeficitAbstraction * getValDivZero(\ - remainingCellAvlWater, - pcr.areatotal(remainingCellAvlWater, allocation_zones), - smallNumber) - # + cellAbstraction += zoneDeficitAbstraction * getValDivZero( + remainingCellAvlWater, + pcr.areatotal(remainingCellAvlWater, allocation_zones), + smallNumber, + ) + # # extraAllocation to minimize numerical errors: - zoneDeficitAllocation = pcr.max(0.0,\ - pcr.areatotal(cellAbstraction, allocation_zones) -\ - pcr.areatotal(cellAllocation , allocation_zones)) + zoneDeficitAllocation = pcr.max( + 0.0, + pcr.areatotal(cellAbstraction, allocation_zones) + - 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, 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 (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: @@ -1612,43 +2117,49 @@ filecache[ncFile] = f # last datetime - last_datetime_year = findLastYearInNCTime(f.variables['time']) - + last_datetime_year = findLastYearInNCTime(f.variables["time"]) + return last_datetime_year - + + def findLastYearInNCTime(ncTimeVariable): # last datetime - last_datetime = nc.num2date(ncTimeVariable[len(ncTimeVariable) - 1],\ - ncTimeVariable.units,\ - ncTimeVariable.calendar) - + last_datetime = nc.num2date( + ncTimeVariable[len(ncTimeVariable) - 1], + ncTimeVariable.units, + ncTimeVariable.calendar, + ) + return last_datetime.year + def findFirstYearInNCTime(ncTimeVariable): # first datetime - first_datetime = nc.num2date(ncTimeVariable[0],\ - ncTimeVariable.units,\ - ncTimeVariable.calendar) - + first_datetime = nc.num2date( + ncTimeVariable[0], ncTimeVariable.units, ncTimeVariable.calendar + ) + return first_datetime.year -def cmd_line(command_line,using_subprocess = True): - msg = "Call: "+str(command_line) +def cmd_line(command_line, using_subprocess=True): + + msg = "Call: " + str(command_line) logger.debug(msg) - + co = command_line if using_subprocess: - cOut,err = subprocess.Popen(co, stdout=subprocess.PIPE,stderr=open('/dev/null'),shell=True).communicate() + cOut, err = subprocess.Popen( + co, stdout=subprocess.PIPE, stderr=open("/dev/null"), shell=True + ).communicate() else: os.system(co) -def plot_variable(pcr_variable, filename = "test.map"): +def plot_variable(pcr_variable, filename="test.map"): + pcr.report(pcr_variable, filename) - cmd = 'aguila '+str(filename) + cmd = "aguila " + str(filename) os.system(cmd) - -