Please update the following lines using xarray slicing:
elif areatype.lower()=="point":
catmask=np.zeros((rainprop.subdimensions))
yind=np.where((latrange-ptlat)>0)[0][-1]
xind=np.where((ptlon-lonrange)>0)[0][-1]
if xind==0 or yind==0:
sys.exit('the point you defined is too close to the edge of the box you defined!')
else:
catmask[yind,xind]=1.0
The point is that xarray should identify the grid cell that is closest to the lat/lon defined by ptlat and ptlon.
Likewise, please replace the following using xarray slicing:
elif areatype.lower()=="box":
finelat=np.arange(latrange[0],latrange[-1]-rainprop.spatialres[1]+rainprop.spatialres[0]/1000,-rainprop.spatialres[1]/25)
finelon=np.arange(lonrange[0],lonrange[-1]+rainprop.spatialres[0]-rainprop.spatialres[0]/1000,rainprop.spatialres[0]/25)
subindy=np.logical_and(finelat>boxarea[2]+rainprop.spatialres[1]/1000,finelat<boxarea[3]+rainprop.spatialres[1]/1000)
subindx=np.logical_and(finelon>boxarea[0]-rainprop.spatialres[0]/1000,finelon<boxarea[1]-rainprop.spatialres[0]/1000)
tx,ty=np.meshgrid(subindx,subindy)
catmask=np.array(np.logical_and(tx==True,ty==True),dtype='float32')
if len(finelat[subindy])<25 and len(finelon[subindx])<25:
print('WARNING: you set POINTAREA to "box", but the box is smaller than a single pixel. This is not advised. Either set POINTAREA to "point" or increase the size of the box.')
if len(finelat[subindy])==1 and len(finelon[subindx])==1:
catmask=np.zeros((rainprop.subdimensions))
yind=np.where((latrange-ptlat)>0)[0][-1]
xind=np.where((ptlon-lonrange)>0)[0][-1]
if xind==0 or yind==0:
sys.exit('the point you defined is too close to the edge of the box you defined!')
else:
catmask[yind,xind]=1.0
else:
def block_mean(ar, fact):
assert isinstance(fact, int), type(fact)
sx, sy = ar.shape
X, Y = np.ogrid[0:sx, 0:sy]
regions = sy//fact * (X//fact) + Y//fact
res = ndimage.mean(ar, labels=regions, index=np.arange(regions.max() + 1))
res.shape = (sx//fact, sy//fact)
return res
catmask=block_mean(catmask,25) # this scheme is a bit of a numerical approximation but I doubt it makes much practical difference
else:
sys.exit("unrecognized area type!")
This should identify the "box" defined by 'boxarea'.
Please update the following lines using xarray slicing: elif areatype.lower()=="point": catmask=np.zeros((rainprop.subdimensions)) yind=np.where((latrange-ptlat)>0)[0][-1] xind=np.where((ptlon-lonrange)>0)[0][-1]
if xind==0 or yind==0: sys.exit('the point you defined is too close to the edge of the box you defined!') else: catmask[yind,xind]=1.0 The point is that xarray should identify the grid cell that is closest to the lat/lon defined by ptlat and ptlon.
Likewise, please replace the following using xarray slicing: elif areatype.lower()=="box": finelat=np.arange(latrange[0],latrange[-1]-rainprop.spatialres[1]+rainprop.spatialres[0]/1000,-rainprop.spatialres[1]/25) finelon=np.arange(lonrange[0],lonrange[-1]+rainprop.spatialres[0]-rainprop.spatialres[0]/1000,rainprop.spatialres[0]/25)
This should identify the "box" defined by 'boxarea'.