Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bugfix for regridding where coarse and fine resolution bounds do not match #31

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
67 changes: 45 additions & 22 deletions model/virtualOS.py
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,7 @@ def singleTryNetcdf2PCRobjCloneWithoutTime(ncFile, varName,\
# 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])
xIdxSta = int(np.where(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0][0])

#~ xIdxSta = int(np.where(np.abs(f.variables['lon'][:] - (xULClone - cellsizeInput/2)) == minX)[0][0])
#~ # see: https://github.com/UU-Hydro/PCR-GLOBWB_model/pull/13
Expand All @@ -203,14 +203,22 @@ def singleTryNetcdf2PCRobjCloneWithoutTime(ncFile, varName,\

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])
yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0][0])

#~ yIdxSta = int(np.where(np.abs(f.variables['lat'][:] - (yULClone - cellsizeInput/2)) == minY)[0][0])
#~ # see: https://github.com/UU-Hydro/PCR-GLOBWB_model/pull/13

#~ yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(factor)))

xULCrop = f.variables['lon'][xIdxSta]-0.5*cellsizeInput
yULCrop = f.variables['lat'][yIdxSta]+0.5*cellsizeInput
if abs(xULClone - xULCrop) > cellsizeClone * 1e-2 or abs(yULClone - yULCrop) > cellsizeClone * 1e-2:
logger.error(f'The the cropped input data corner (XUIL: {xULCrop} - yUL: {yULCrop} - cellsize {cellsizeInput}) does not align with the clone map (XUIL: {xULClone} - yUL: {yULClone} - cellsize {cellsizeClone}). ' +
'This usually occurs if the clone-map does not align with the input data. ' +
'Please make sure to have the clone-map corners align with whole degrees of latitude and longitude.')


cropData = f.variables[varName][yIdxSta:yIdxEnd,xIdxSta:xIdxEnd]

if factor > 1: logger.debug('Resample: input cell size = '+str(float(cellsizeInput))+' ; output/clone cell size = '+str(float(cellsizeClone)))
Expand All @@ -228,18 +236,22 @@ def singleTryNetcdf2PCRobjCloneWithoutTime(ncFile, varName,\

# convert to PCR object and close f
if specificFillValue != None:
outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData2FinerGrid(factor, cropData, float(specificFillValue)), \
float(specificFillValue))
fillValue = float(specificFillValue)
else:
try:
outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData2FinerGrid(factor, cropData, float(f.variables[varName]._FillValue)), \
float(f.variables[varName]._FillValue))
fillValue = float(f.variables[varName]._FillValue)
except:
outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData2FinerGrid(factor, cropData, float(f.variables[varName].missing_value)), \
float(f.variables[varName].missing_value))
fillValue = float(f.variables[varName].missing_value)

regridData = regridData2FinerGrid(factor, cropData, fillValue)
if regridData.shape[0] % factor != 0 or regridData.shape[1] % factor != 0:
logger.error(f'The the regridded input data shape ({regridData.shape}) does not align with the clone map shape (rows: {rowsClone} - cols: {colsClone}). ' +
'This usually occurs if the clone-map does not align with the input data. ' +
'Please make sure to have the clone-map corners align with whole degrees of latitude and longitude.')

outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData, \
fillValue)

#~ # debug:
#~ pcr.report(outPCR,"tmp.map")
Expand Down Expand Up @@ -804,7 +816,7 @@ def singleTryNetcdf2PCRobjClone(ncFile,\
# 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])
xIdxSta = int(np.argmin(abs(f.variables['lon'][:] - (xULClone + 0.5*cellsizeInput)) == minX)[0][0])

#~ xIdxSta = int(np.where(np.abs(f.variables['lon'][:] - (xULClone - cellsizeInput/2)) == minX)[0][0])
#~ # see: https://github.com/UU-Hydro/PCR-GLOBWB_model/pull/13
Expand All @@ -814,14 +826,21 @@ def singleTryNetcdf2PCRobjClone(ncFile,\

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])
yIdxSta = int(np.where(abs(f.variables['lat'][:] - (yULClone - 0.5*cellsizeInput)) == minY)[0][0])

#~ yIdxSta = int(np.where(np.abs(f.variables['lat'][:] - (yULClone - cellsizeInput/2)) == minY)[0][0])
#~ # see: https://github.com/UU-Hydro/PCR-GLOBWB_model/pull/13

#~ yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(cellsizeInput/cellsizeClone)))
yIdxEnd = int(math.ceil(yIdxSta + rowsClone /(factor)))

xULCrop = f.variables['lon'][xIdxSta]-0.5*cellsizeInput
yULCrop = f.variables['lat'][yIdxSta]+0.5*cellsizeInput
if abs(xULClone - xULCrop) > cellsizeClone * 1e-2 or abs(yULClone - yULCrop) > cellsizeClone * 1e-2:
logger.error(f'The the cropped input data corner (XUIL: {xULCrop} - yUL: {yULCrop} - cellsize {cellsizeInput}) does not align with the clone map (XUIL: {xULClone} - yUL: {yULClone} - cellsize {cellsizeClone}). ' +
'This usually occurs if the clone-map does not align with the input data. ' +
'Please make sure to have the clone-map corners align with whole degrees of latitude and longitude.')

# retrieve data from netCDF for slice

if f.variables[varName].ndim == 4:
Expand Down Expand Up @@ -851,18 +870,22 @@ def singleTryNetcdf2PCRobjClone(ncFile,\

# convert to PCR object and close f
if specificFillValue != None:
outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData2FinerGrid(factor, cropData, float(specificFillValue)), \
float(specificFillValue))
fillValue = float(specificFillValue)
else:
try:
outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData2FinerGrid(factor, cropData, float(f.variables[varName]._FillValue)), \
float(f.variables[varName]._FillValue))
fillValue = float(f.variables[varName]._FillValue)
except:
outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData2FinerGrid(factor, cropData, float(f.variables[varName].missing_value)), \
float(f.variables[varName].missing_value))
fillValue = float(f.variables[varName].missing_value)

regridData = regridData2FinerGrid(factor, cropData, fillValue)
if regridData.shape[0] % factor != 0 or regridData.shape[1] % factor != 0:
logger.error(f'The the regridded input data shape ({regridData.shape}) does not align with the clone map shape (rows: {rowsClone} - cols: {colsClone}). ' +
'This usually occurs if the clone-map does not align with the input data. ' +
'Please make sure to have the clone-map corners align with whole degrees of latitude and longitude.')

outPCR = pcr.numpy2pcr(pcr.Scalar, \
regridData, \
fillValue)

#~ pcr.aguila(outPCR)

Expand Down