def region(xds, name='REGION1', ra=None, dec=None, pixels=None, pol=-1,
    channels=-1):
    """
    Create a new region Data variable in the Dataset 

    
    .. note:: This function currently only supports rectangles and integer pixel boundaries
    
    Parameters
    ----------
    xds : xarray.core.dataset.Dataset
        input image dataset
    name : str
        dataset variable name for region, overwrites if already present
    ra : list
        right ascension coordinate range in the form of [min, max]. Default None means all
    dec : list
        declination coordinate range in the form of [min, max]. Default None means all
    pixels : array_like
        array of shape (N,2) containing pixel box. OR'd with ra/dec
    pol : int or list
        polarization dimension(s) to include in region.  Default of -1 means all
    channels : int or list
        channel dimension(s) to include in region.  Default of -1 means all
        
    Returns
    -------
    xarray.core.dataset.Dataset
        New Dataset
    """
    import numpy as np
    import xarray as xr
    if not name.strip():
        name = 'REGIONX'
    if ra is None:
        ra = [0.0, 0.0]
    if dec is None:
        dec = [0.0, 0.0]
    if pixels is None:
        pixels = np.zeros((1, 2), dtype=int) - 1
    pixels = np.array(pixels, dtype=int)
    if pixels.ndim != 2 or pixels.shape[1] != 2:
        print('ERROR: pixels parameter not a (N,2) array')
        return None
    pol = np.array(np.atleast_1d(pol), dtype=int)
    if pol[0] == -1:
        pol = list(range(len(xds['pol'])))
    channels = np.array(np.atleast_1d(channels), dtype=int)
    if channels[0] == -1:
        channels = list(range(len(xds['chan'])))
    region = xr.ones_like(xds.IMAGE, dtype=bool).where((xds.right_ascension >
        np.min(ra)) & (xds.right_ascension < np.max(ra)) & (xds.declination >
        np.min(dec)) & (xds.declination < np.max(dec)), False)
    region = region | xr.ones_like(xds.IMAGE, dtype=bool).where((xds.d0 >
        np.min(pixels[:, (0)])) & (xds.d0 < np.max(pixels[:, (0)])) & (xds.
        d1 > np.min(pixels[:, (1)])) & (xds.d1 < np.max(pixels[:, (1)])), False
        )
    region = region.where(xds.pol.isin(xds.pol[pol]), False)
    region = region.where(xds.chan.isin(xds.chan[channels]), False)
    xds = xds.assign(dict([(name, region)]))
    return xds
