Bin data on a regular latitude and longitude grid

bin_lat_lon.py

def bin_lat_lon(data, lon_min=0., lon_max=360., lon_res=1., 
                lat_min=-80., lat_max=80, lat_res=1.):
    """
    Bin unidimensional data onto a regular latidude and longitude grid
    by computing the median value and the number of corresponding 
    observations
    
    Paraneters
    ----------
    data : xarray.DataArray
        The data stored in a DataArray object
    lon_min : float, optional
        The western longitudinal boundary (default is 0.)
    lon_max : float, optional
        The eastern longitudinal boundary (default is 360.)
    lon_res : float, optional
        The longitudinal resolution of the grid (default is 1.)
    lat_min : float, optional
        The southern latitudinal boundary (default is -80.)
    lat_max : float, optional
        The northern latitudinal boundary (default is 80.)
    lat_res : float, optional
        The latitudinal resolution of the grid (default is 1.)
    
    Returns
    -------
    res : xarray.Dataset
        A gridded dataset containing the median value and the number of 
        observation per bins    
    """
    lon_bins = np.arange(lon_min, lon_max, lon_res)
    lon_labels = lon_bins[:-1] - np.diff(lon_bins) / 2
    lat_bins = np.arange(lat_min, lat_max, lat_res)
    lat_labels = lat_bins[:-1] - np.diff(lat_bins) / 2
    mean_values = []
    total_nobs = []
    lat_values = []
    # Loop over latitudes
    for i, ds in list(data.groupby_bins('nav_lat', lat_bins, 
                                        labels=lat_labels, 
                                        include_lowest=True)):
        try:
            group = ds.groupby_bins('nav_lon', lon_bins, 
                                    labels=lon_labels, 
                                    include_lowest=True)
            bins = group.median().sortby('lon_bins')
            nobs = group.count().sortby('lon_bins')
            mean_values.append(bins)
            total_nobs.append(nobs)
            lat_values.append(i)
        except (ValueError, StopIteration):
            dummy_array = xr.DataArray(np.full(len(lon_labels), np.nan), 
                                       dims='lon_bins', 
                                       coords={'lon_bins': 
                                               ('lon_bins', lon_labels)
                                              }
                                      )
            mean_values.append(dummy_array)
            lat_values.append(i)
    res_bins = (xr.concat(mean_values, dim='nav_lat')
                  .assign_coords(lat=lat_values)
                  .rename({'lon_bins': 'nav_lon'})
                  .sortby('nav_lat')
               )
    res_obs = (xr.concat(total_nobs, dim='nav_lat')
                 .assign_coords(lat=lat_values)
                 .rename({'lon_bins': 'nav_lon'})
                 .sortby('nav_lat')
              )
    res = xr.Dataset({data.name:res_bins, 'nobs':res_obs})
    return res