diff --git a/src/pyobs/csBin.py b/src/pyobs/csBin.py index 3acbcb9..bbf48d1 100644 --- a/src/pyobs/csBin.py +++ b/src/pyobs/csBin.py @@ -221,7 +221,7 @@ def setIndices (self, lons, lats): return - def binObs ( self, obs ): + def binObs ( self, obs, average = True ): """ Given a list of longitude and latitudes in (lons,lats), bin list of observations *obs* on the cubed sphere, @@ -257,8 +257,12 @@ def binObs ( self, obs ): # --------- gObs = MAPL_UNDEF * np.ones(face_size, dtype=float) K = nObs>0 - gObs[K] = aObs[K] / nObs[K] - + + if (average) : + gObs[K] = aObs[K] / nObs[K] + else : + gObs[K] = aObs[K] + csObs[f] = np.reshape(gObs, (self.Xdim, self.Xdim)) # All done # -------- diff --git a/src/pyobs/llBin.py b/src/pyobs/llBin.py new file mode 100644 index 0000000..6145555 --- /dev/null +++ b/src/pyobs/llBin.py @@ -0,0 +1,183 @@ +""" + + Bin observations on LatLon grid. + +""" + +import numpy as np +import bisect +import xarray as xr + +from .constants import MAPL_UNDEF + +class llBinError(Exception): + """ + Defines general exception errors. + """ + def __init__(self, value): + self.value = value + def __str__(self): + return repr(self.value) + +class LLBIN(object): + + def __init__(self, **kwargs): + """ + Supporting keywords: GridName,lon_bnds, lat_bnds. + GridName: It should have the standard format : POLEiiiXjjj-DATELINE + POLE : PE or PC + iii : number of lon + x : the separator + jjj : number of lat + DATELINE : DE or DC + + For example, PE360x180-DE represents a lat-lon grid with 360 lons, 180 lats, pole on edge and Dateline on edge + + lon_bnds: longitude bounds, an ordered array from -180. to 180. + lat_bnds: latitude bounds, an ordered array from -90. to 90. + """ + + kw_ = {k.upper():v for k,v in kwargs.items()} + + if 'GRIDNAME' in kw_.keys(): + grid = kw_['GRIDNAME'].upper() + pole = grid[:2] + x = grid.find('X') + dateline = grid[-2:] + nlon = int(grid[2:x]) + nlat = int(grid[x+1:-3]) + self.pole= pole + self.dateline=dateline + self.Xdim = nlon + self.Ydim = nlat + elif ( 'LAT_BNDS' in kw_.keys() and 'LON_BNDS' in kw_.keys()): + self.lon_bnds = lon_bnds + self.lat_bnds = lat_bnds + self.pole='unknown' + self.dateline='unknown' + self.Xdim = len(self.lon_bnds)-1 + self.Ydim = len(self.lat_bnds)-1 + else: + raise llBinError('Init keywork is not correct') + + def setIndices (self, lons, lats): + """ + Given a list of longitude and latitudes in degree (lons,lats), + find the cell indices. + """ + I = np.zeros(len(lons)) + J = np.zeros(len(lats)) + # get I indices + if (self.dateline =='DE'): + dlon = 360./self.Xdim + I = np.floor((lons+180.)/dlon).astype(int) + elif (self.dateline =='DC'): + dlon = 360./self.Xdim + lons = lons + dlon/2. + lons = np.where(lons >= 180.0, lons - 360., lons) + I = np.floor((lons+180.0)/dlon).astype(int) + elif (self.dateline =='unknown'): + I = np.array([ bisect.bisect(self.lon_bnds, lon) for lon in lons]) - 1 + # get I indices + if (self.pole == 'PE'): + dlat = 180./self.Ydim + J = np.floor((lats + 90.)/dlat).astype(int) + elif (self.pole == 'PC'): + dlat = 180./(self.Ydim-1) + J = np.floor((lats + 90.+ dlat/2.)/dlat).astype(int) + elif (self.pole == 'unknown'): + J = np.array([ bisect.bisect(self.lat_bnds, lat) for lat in lats]) - 1 + # m + I = np.where( I == self.Xdim, 0, I) + J = np.where( J == self.Ydim, J-1, J) + + self.IJ = I + J*self.Xdim + + return + + def binObs ( self, obs, average = True ): + """ + Given a list of longitude and latitudes in (lons,lats), + bin list of observations *obs* on the lat-lon grid, + returning gridded observations *llObs*. + """ + + llObs = np.zeros((self.Xdim, self.Ydim), dtype=float) + grid_size = self.Xdim * self.Ydim + aObs = np.zeros(grid_size, dtype=float) + nObs = np.zeros(grid_size, dtype=int) + + if (len(obs.shape) > 1): + obs = np.ravel(obs) + + # Accumulator, counter + # -------------------- + + # Accumulate observations lying on this grid + # ------------------------------------------------------ + try: + IJ = self.IJ + except: + raise llBinError('Indices (F,IJ) not yet set.') + + np.add.at(nObs, IJ, 1) + np.add.at(aObs, IJ, obs) + + # Normalize + # --------- + gObs = MAPL_UNDEF * np.ones(grid_size, dtype=float) + K = nObs>0 + + if (average) : + gObs[K] = aObs[K] / nObs[K] + else : + gObs[K] = aObs[K] + + llObs = np.reshape(gObs, (self.Xdim, self.Xdim)) + # All done + # -------- + return llObs + +#........................................................................... +if __name__ == '__main__': + + llbin = LLBIN(GridName='PE360x180-DE') + nlon = 360 + nlat = 180 + dlat = 180./nlat + dlon = 360./nlon + + lats = xr.DataArray([ -90.]) + lons = xr.DataArray([ -180.]) + llbin.setIndices(lons, lats) + print(llbin.IJ[0] == 0) + + lats = xr.DataArray([ -90. + 1.5*dlat ]) + lons = xr.DataArray([ -180.]) + llbin.setIndices(lons, lats) + print(llbin.IJ[0] == 360) + + lats = xr.DataArray([ -90. +1.5*dlat ]) + lons = xr.DataArray([ -180. +1.5*dlon ]) + llbin.setIndices(lons, lats) + print(llbin.IJ[0] == 361) + + lats = xr.DataArray([ -90. +10.5*dlat ]) + lons = xr.DataArray([ -180. +10.5*dlon ]) + llbin.setIndices(lons, lats) + print(llbin.IJ[0] == 3610) + + lon_bnds = xr.DataArray([ -180.0+i*dlon for i in range(nlon+1)]) + lat_bnds = xr.DataArray([ -90.0 +i*dlat for i in range(nlat+1)]) + + llbin = LLBIN(lon_bnds=lon_bnds, lat_bnds=lat_bnds) + + lats = xr.DataArray([ -90. + 0.1*dlat]) + lons = xr.DataArray([ -180. + 0.1*dlon]) + llbin.setIndices(lons, lats) + print(llbin.IJ[0] == 0 ) + + lats = xr.DataArray([ -90. +10.5*dlat ]) + lons = xr.DataArray([ -180. +10.5*dlon ]) + llbin.setIndices(lons, lats) + print(llbin.IJ[0] == 3610)