Repository: ks905383/xagg

Issue:
0.3.3: `numba` implementation doesn't do secondary weights yet!
Self-explanatory. For some reason, this wasn't caught by the unit tests. 

TO-DO: 
- Make sure unit tests correctly break when additional weights aren't included
- Add additional weights support to `numba` implementation

Discussion:

Solution patch:
diff --git a/ci/environment-py3.10.yml b/ci/environment-py3.10.yml
index c6b1d96..77b5bf6 100644
--- a/ci/environment-py3.10.yml
+++ b/ci/environment-py3.10.yml
@@ -16,6 +16,7 @@ dependencies:
   - esmf >= 8.1.0
   - esmpy >= 8.1.0 
   - pytables
+  - dask
   ##############
   - pytest
   - pip:
diff --git a/ci/environment-py3.11.yml b/ci/environment-py3.11.yml
index 252817c..90351d4 100644
--- a/ci/environment-py3.11.yml
+++ b/ci/environment-py3.11.yml
@@ -16,6 +16,7 @@ dependencies:
   - esmf >= 8.1.0
   - esmpy >= 8.1.0 
   - pytables
+  - dask
   ##############
   - pytest
   - pip:
diff --git a/ci/environment-py3.12.numba.yml b/ci/environment-py3.12.numba.yml
index e5bf22f..38bca61 100644
--- a/ci/environment-py3.12.numba.yml
+++ b/ci/environment-py3.12.numba.yml
@@ -17,6 +17,7 @@ dependencies:
   - esmpy >= 8.1.0 
   - pytables
   - numba
+  - dask 
   ##############
   - pytest
   - pip:
diff --git a/ci/environment-py3.12.plotfuncs.yml b/ci/environment-py3.12.plotfuncs.yml
index 8270ee8..28b3297 100644
--- a/ci/environment-py3.12.plotfuncs.yml
+++ b/ci/environment-py3.12.plotfuncs.yml
@@ -19,6 +19,7 @@ dependencies:
   - cartopy
   - matplotlib
   - cmocean
+  - dask
   ##############
   - pytest
   - pip:
diff --git a/ci/environment-py3.12.yml b/ci/environment-py3.12.yml
index 94c272d..8acb55e 100644
--- a/ci/environment-py3.12.yml
+++ b/ci/environment-py3.12.yml
@@ -16,6 +16,7 @@ dependencies:
   - esmf >= 8.1.0
   - esmpy >= 8.1.0 
   - pytables
+  - dask
   ##############
   - pytest
   - pip:
diff --git a/ci/environment-py3.9.yml b/ci/environment-py3.9.yml
index 1968ed8..4efd7ec 100644
--- a/ci/environment-py3.9.yml
+++ b/ci/environment-py3.9.yml
@@ -16,6 +16,7 @@ dependencies:
   - esmf >= 8.1.0
   - esmpy >= 8.1.0 
   - pytables
+  - dask
   ##############
   - pytest
   - pip:
diff --git a/xagg/auxfuncs.py b/xagg/auxfuncs.py
index 89f86d6..86d175e 100644
--- a/xagg/auxfuncs.py
+++ b/xagg/auxfuncs.py
@@ -70,6 +70,133 @@ def list_or_first(ser):
     else:
         return lis
 
+def _diagnose_nans(da,other_dims,
+                   _return_somelocs = False):
+    ''' Flag if partial nans in other dims
+
+    Parameters
+    ------------
+    da : :py:class:`xr.DataArray`
+
+    other_dims : :py:class:`list` of :py:meth:`str`
+        Dimensions over which to test nan prevalence
+
+    Returns
+    ------------
+    nanflags : :py:class:`dict`
+        `nanflags['all']` : A dictionary with keys of dimension names 
+                            and values True if any of the coordinates
+                            in that dimension has all nans in the `da`
+        `nanflags['some']` : A dictionary with keys of dimension names
+                             and values True if some coordinates in 
+                             that dimension have nans that other coordinates
+                             do not
+
+
+    '''
+    # Copy, to make sure the original dataarray is untouched
+    da = da.copy()
+
+    # Flag if all values of the variable are nan for some specific coordinate
+    # along a dimension (this is not per se a problem, but good to have a 
+    # flag if needed) 
+    allnan_flags = {dim:da.isnull().all(dim = [d for d in da.sizes if d != dim]).any().values
+                    for dim in other_dims}
+
+    # Flag if some values of the variable are nan for some specific coordinate
+    # along a dimension (this _could_ be a problem, because it means that 
+    # the weights will be different between different values of that coordinate
+    # First, drop all-nan slices, since they're captured above and not as much 
+    # of a problem
+    for dim in other_dims: 
+        da = da.dropna(dim,how='all')
+    # Now, find for which `locs` there are inconsistent nans 
+    somenan_flags = {dim:(da.isnull().any(dim = dim) != 
+                          da.isnull().all(dim = dim)).any().values
+                         if da.sizes[dim] != 0 else False
+                     for dim in other_dims}
+
+    # Combine into single dict
+    nanflags = {'all':allnan_flags,
+                'some':somenan_flags}
+
+    if _return_somelocs:
+        somelocs = {dim:(da.isnull().any(dim = dim) != 
+                          da.isnull().all(dim = dim)).any(dim = [d for d in other_dims if d != dim])
+                     for dim in other_dims}
+
+        return nanflags,somelocs
+    else:
+        return nanflags
+ 
+class SomeNanWarning(UserWarning):
+    ''' Warns when pixels are inconsistently nan along a dimension 
+    '''
+    pass
+
+def _warn_ifsomenans(ds,var,dims,
+                    _warn_trigger_partialnan=True,
+                    _return_somelocs=False):
+    ''' Wrapper for SomeNanWarning with standard text.
+    Triggers warning if `_warn_trigger_partialnan` is True and 
+    `_diagnose_nans(ds[var],other_dims)` returns True for at least
+    one dim for `nanflags['some']` 
+
+    ## TODO - is there a way to flag which poly_idxs this would 
+    affect? Should be, since this is used in xa.aggregate, that has
+    that info. Basically, can do a [is in] across pix_idxs to flag
+    poly_idxs that should optionally be nan in this case 
+    '''
+    if type(ds) == xr.core.dataset.Dataset:
+        ds = ds[var].copy()
+    
+    if _return_somelocs:
+        raise NotImplementedError
+
+        nanflags,somelocs = _diagnose_nans(ds,dims,
+             _return_somelocs=True)
+
+    else:
+        nanflags = _diagnose_nans(ds,dims)
+
+
+    if (_warn_trigger_partialnan and 
+        np.any([flag for dim,flag in nanflags['some'].items()])):
+        dims_with_somenans = [dim for dim,flag in nanflags['some'].items() if flag]
+        
+        warnings.warn('One or more grid cells in variable '+var+' have inconsistent nans along the dimension(s) '+
+                       ', '.join(dims_with_somenans)+' (i.e., one or more grid cells are nan for some but not all coordinates of the dimension(s))'+
+                       '. This means that grid cell weights will be different for different coordinates along '+
+                       ', '.join(dims_with_somenans)+'. The aggregation '+
+                       'calculation may therefore be incorrect, since it aggregates over different grid cells for '+
+                       ' the same polgyon for different coordinates of the dimension(s)'+', '.join(dims_with_somenans)+'.',
+                       SomeNanWarning)
+
+        _warn_trigger_partialnan = False
+
+    if _return_somelocs:
+        from functools import reduce
+        import operator
+
+        # Get all locations that are at partially nan in one direction
+        bad_locs = reduce(operator.or_, somelocs.values())
+        bad_locs = bad_locs.where(bad_locs,drop=True)
+
+        # Get their lat, lon locations
+        bad_coords = [(lat,lon) for lat,lon in zip(bad_locs.lat.values,bad_locs.lon.values)]
+        # Turn into set (for faster indexing)
+        bad_coords = set(bad_coords)
+
+        # Check each list of polygon grid cell overlap coordinates if they contain any bad coordinate
+        # This will need to go into aggregate()
+        #poly_badloc_mask = wm.agg.coords.apply(lambda x: bool(bad_coords.intersection(x)))
+        # and then a later, something that replaces the output in those locations with nan
+
+        return _warn_trigger_partialnan,bad_coords
+    else:
+        return _warn_trigger_partialnan
+   
+
 
 def fix_ds(ds,var_cipher = {'latitude':{'latitude':'lat','longitude':'lon'},
                             'Latitude':{'Latitude':'lat','Longitude':'lon'},
diff --git a/xagg/core.py b/xagg/core.py
index 1f61749..948150b 100644
--- a/xagg/core.py
+++ b/xagg/core.py
@@ -18,7 +18,7 @@ try:
 except ImportError:
     _has_numba = False
 
-from . auxfuncs import (find_rel_area,normalize,fix_ds,get_bnds,subset_find,list_or_first)
+from . auxfuncs import (find_rel_area,normalize,fix_ds,get_bnds,subset_find,list_or_first,_warn_ifsomenans)
 from . classes import (weightmap,aggregated)
 from . options import get_options
 
@@ -570,30 +570,35 @@ if _has_numba:
         return vals
 
     @njit
-    def numba_aggregate(data,idxs,weights):
+    def numba_aggregate(data,idxs,rel_area,add_weights=None):
+        if add_weights is None:
+            add_weights = np.ones(len(data))
+        
         # Since idxs, rel_area are stored in a dataarray,
         # in the weightmap, they're padded with nans. 
         # Remove those. 
         idxs = idxs[~np.isnan(idxs)]
-        weights = weights[~np.isnan(weights)]
-
+        rel_area = rel_area[~np.isnan(rel_area)]
+    
         # Get the indices of the data needed
         data = numba_indexer(data,idxs)
+        add_weights = numba_indexer(add_weights,idxs)
 
+        # Normalize area and additional weights to a single
+        # weights vector
+        weights = rel_area * add_weights
+    
         # Make sure to do aggregation only on non-nan pixels
         nans = np.isnan(data)
         data = data[~nans]
         weights = weights[~nans]
-
+    
         if len(data) == 0:
             agg = np.nan
         else:
             agg = np.sum(data*weights)/np.sum(weights)
-
+    
         return agg
-
-    def wrapper_numba_aggregate(data,idxs,rel_area):
-        return agggregate_numba(data,idxs,rel_area)
 #--------------------------------------------
 
 def aggregate(ds,wm,impl=None,silent=None):
@@ -614,8 +619,8 @@ def aggregate(ds,wm,impl=None,silent=None):
     
     Parameters
     ---------------
-    ds : :class:`xarray.Dataset`
-        an :class:`xarray.Dataset` containing one or more
+    ds : :py:class:`xarray.Dataset`
+        an :py:class:`xarray.Dataset` containing one or more
         variables with dimensions lat, lon (and possibly
         more). The dataset's geographic grid has to 
         include the lat/lon coordinates used in 
@@ -623,9 +628,9 @@ def aggregate(ds,wm,impl=None,silent=None):
         :func:`xagg.core.get_pixel_overlaps` (and saved in 
         ``wm['source_grid']``)
                
-    wm : :class:`xagg.classes.weightmap`
+    wm : :py:class:`xagg.classes.weightmap`
         the output to :func:`xagg.core.get_pixel_overlaps`; a 
-        :class:`xagg.classes.weightmap` object containing 
+        :py:class:`xagg.classes.weightmap` object containing 
 
         - ``['agg']``
             a dataframe, 
@@ -638,7 +643,7 @@ def aggregate(ds,wm,impl=None,silent=None):
             were calculated (and on which the linear indices 
             are based)
 
-    impl : :class:str (def: ``'for_loop'``) (set by :py:meth:`xa.set_options`)
+    impl : :py:class:str (def: ``'for_loop'``) (set by :py:meth:`xa.set_options`)
         which aggregation calculation method to use, either of: 
 
         - ``'for_loop'`` 
@@ -725,18 +730,10 @@ def aggregate(ds,wm,impl=None,silent=None):
 
                 # Check first if any nans are "complete" (meaning that a pixel 
                 # either has values for each step, or nans for each step - if
-                # there are random nans within a pixel, throw a warning)
-                if (_warn_trigger_partialnan and 
-                    (not xr.Dataset.equals(np.isnan(var_array).any(other_dims),
-                                           np.isnan(var_array).all(other_dims)))):
-                    warnings.warn('One or more pixels in variable '+var+' have *some* nans in the dimension(s) '+
-                                           ', '.join(other_dims)+
-                                           '. The code can currently only deal with pixels for which the '+
-                                           'pixel is *entirely* nan in all dimensions (or has no nans), however there is currently no '+
-                                           ' support for data in which pixels have only some nan values. The aggregation '+
-                                           'calculation is likely incorrect, since the weights will be different across different '+
-                                           'coordinates in the dimension(s) '+', '.join(other_dims)+'.')
-                    _warn_trigger_partialnan = False
+                # there are random nans along non-location dimensions for the 
+                # same grid cell, throw a warning)
+                _warn_trigger_partialnan=_warn_ifsomenans(var_array,var,other_dims,
+                                                            _warn_trigger_partialnan)
 
 
                 # multiply percent-overlaps by user-supplied weights 
@@ -765,6 +762,8 @@ def aggregate(ds,wm,impl=None,silent=None):
                 data_dict[var] = aggregated_array
 
         ds_combined = xr.Dataset(data_dict)  
+        print('ds_combined:')
+        print(ds_combined)
 
         for var in ds:
             if ('bnds' not in ds[var].sizes) and ('loc' in ds[var].sizes):
@@ -799,18 +798,11 @@ def aggregate(ds,wm,impl=None,silent=None):
                         other_dims = [k for k in np.atleast_1d(ds[var].sizes) if k != 'loc']
                         # Check first if any nans are "complete" (meaning that a pixel 
                         # either has values for each step, or nans for each step - if
-                        # there are random nans within a pixel, throw a warning)
-                        if (_warn_trigger_partialnan and 
-                            (not xr.Dataset.equals(np.isnan(ds[var].isel(loc=wm.agg.iloc[poly_idx,:].pix_idxs)).any(other_dims),
-                                                   np.isnan(ds[var].isel(loc=wm.agg.iloc[poly_idx,:].pix_idxs)).all(other_dims)))):
-                            warnings.warn('One or more pixels in variable '+var+' have *some* nans in the dimension(s) '+
-                                           ', '.join(other_dims)+
-                                           '. The code can currently only deal with pixels for which the '+
-                                           'pixel is *entirely* nan in all dimensions (or has no nans), however there is currently no '+
-                                           ' support for data in which pixels have only some nan values. The aggregation '+
-                                           'calculation is likely incorrect, since the weights will be different across different '+
-                                           'coordinates in the dimension(s) '+', '.join(other_dims)+'.')
-                            _warn_trigger_partialnan = False
+                        # there are random nans along non-location dimensions for the 
+                        # same grid cell, throw a warning)
+                        _warn_trigger_partialnan=_warn_ifsomenans(ds,var,other_dims,
+                                                                  _warn_trigger_partialnan)
+
 
                         if not np.isnan(ds[var].isel(loc=wm.agg.iloc[poly_idx,:].pix_idxs)).all(): 
                             # Get relative areas for the pixels overlapping with this Polygon
@@ -841,26 +833,65 @@ def aggregate(ds,wm,impl=None,silent=None):
         if not _has_numba:
             raise ImportError("impl == 'numba' requires `numba`, which is not installed.")
         # Extract pixel idxs, area overlap for each polygon from weightmap object
-        idxs = xr.concat([xr.Dataset({'idxs':(['idx'],np.atleast_1d(np.array(idxs))),
-                     'rel_area':(['idx'],np.atleast_1d(np.atleast_1d(np.array(areas))[0]))},
-                      coords = {'idx':(['idx'],np.arange(len(np.atleast_1d(idxs))))}) 
-         for idxs,areas in zip(wm.agg.pix_idxs.values,wm.agg.rel_area.values)],
-                  dim=pd.Index(wm.agg.poly_idx.values,name='poly_idx'))
+        pix_idxs = [np.atleast_1d(pix) for pix in wm.agg.pix_idxs.values]
+        rel_area = [np.atleast_1d(area[0]) if area is not np.nan else 
+                    np.atleast_1d(np.nan) for area in wm.agg.rel_area.values]
+        max_pix = max(len(lst) for lst in pix_idxs)  # Max number of pixels per polygon, to pad to
+
+        # Preallocate arrays filled with NaNs
+        n_polys = len(wm.agg)
+        pix_idxs_padded = np.full((n_polys, max_pix), np.nan)
+        rel_area_padded = np.full((n_polys, max_pix), np.nan)
+
+        # Fill the arrays, padding to max number of pixels with nans
+        for i, (pix_idxs, rel_area) in enumerate(zip(pix_idxs, rel_area)):
+            pix_idxs_padded[i, :len(pix_idxs)] = pix_idxs
+            rel_area_padded[i, :len(rel_area)] = rel_area
+
+        # Put into dataset
+        idxs = xr.Dataset({'idxs': (['poly_idx', 'idx'], pix_idxs_padded),
+                          'rel_area': (['poly_idx', 'idx'], rel_area_padded)},
+                          coords={'poly_idx': wm.agg.poly_idx.values,
+                                    'idx': np.arange(max_pix)}
+                        )
+        # Make into dask array 
+        idxs = idxs.chunk({'poly_idx':'auto','idx':-1})
+
+        # Make weights data array, to put into apply_ufunc below
+        weights = xr.DataArray(weights,
+                            dims = ['loc'],
+                            coords = {'loc':ds.loc})
 
         for var in ds:
             if ('bnds' not in ds[var].sizes) and ('loc' in ds[var].sizes):
                 if not silent:
                     print('aggregating '+var+'...')
-                agg = xr.merge([xr.apply_ufunc(numba_aggregate,
+
+                # Get the dimensions of the variable that aren't "loc" (location)
+                other_dims = [k for k in np.atleast_1d(ds[var].sizes) if k != 'loc']
+
+                # Check first if any nans are "complete" (meaning that a pixel 
+                # either has values for each step, or nans for each step - if
+                # there are random nans along non-location dimensions for the 
+                # same grid cell, throw a warning)
+                _warn_trigger_partialnan=_warn_ifsomenans(ds,var,other_dims,
+                                                          _warn_trigger_partialnan)
+
+                # Aggregate
+                agg = xr.apply_ufunc(numba_aggregate,
                                    ds[var],
                                    idxs.idxs,
                                    idxs.rel_area,
-                                   input_core_dims = [['loc'],['idx'],['idx']],
+                                   weights,
+                                   input_core_dims = [['loc'],['idx'],['idx'],['loc']],
                                    output_core_dims = [[]],
                                    vectorize=True,
                                    dask='parallelized',
                                    output_dtypes = [float]).to_dataset(name=var)
-                                for var in ds if (('bnds' not in ds[var].sizes) and ('loc' in ds[var].sizes))])
+
+                # Trigger computation, otherwise it separately triggers computation
+                # for every poly_idx below, with tremendous overhead 
+                agg = agg.compute()
 
                 # Put in `aggregated` format (TODO: should reform this to remove some of the 
                 # nested lists that are here for some reason)


Test patch:
diff --git a/tests/test_auxfuncs.py b/tests/test_auxfuncs.py
index 8ada5dd..13ed7a0 100644
--- a/tests/test_auxfuncs.py
+++ b/tests/test_auxfuncs.py
@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
 import xarray as xr
-from xagg.auxfuncs import (normalize,fix_ds,get_bnds,subset_find)
+from xagg.auxfuncs import (normalize,fix_ds,get_bnds,subset_find,_diagnose_nans)
 
 ##### normalize() tests #####
 def test_normalize():
@@ -310,5 +310,44 @@ def test_subset_find_nomatch():
 		subset_find(ds0,ds1)
 
 
+###### _diagnose_nans tests ######
+dims = ['lat','lon','dim1','dim2']
+da = xr.DataArray(np.random.rand(*([2]*len(dims))),
+                  dims = dims,
+                  coords = {dim:[0,1] for dim in dims})
+
+# Stack, which _diagnose_nans expects (well, not really, but 
+# this is the form that it'll always get it)
+da = da.stack(loc = ['lat','lon'])
+
+def test_diagnose_nans_all(da=da):
+    da = da.copy()
+
+    # Now, make all values along one of the dim1 coordinates nan
+    da.loc[{'dim1':0}] = np.nan
+
+    # Run _diagnose_nans
+    nanflags = _diagnose_nans(da,other_dims = ['dim1','dim2'])
+
+    # Make sure 'all' and only 'all' is triggered, and only for dim 1
+    assert nanflags['all']['dim1']
+    assert not nanflags['all']['dim2']
+    assert np.all(~flag for flag in nanflags['some'])
+
+def test_diagnose_nans_some(da=da):
+    da = da.copy()
+
+    # Now, make one value in a single location alog one of the dim1 
+    # coordinates nan
+    da.loc[{'lon':0,'dim1':0}] = np.nan
+
+    # Run _diagnose_nans
+    nanflags = _diagnose_nans(da,other_dims = ['dim1','dim2'])
+
+    # Make sure 'some' and only 'some' is triggered, and only for dim 1
+    assert nanflags['some']['dim1']
+    assert not nanflags['some']['dim2']
+    assert np.all(~flag for flag in nanflags['all'])
+
 
 
diff --git a/tests/test_core.py b/tests/test_core.py
index 8f5a092..bccac02 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -25,6 +25,7 @@ except ImportError:
     _has_numba = False
 
 from xagg.core import (process_weights,create_raster_polygons,get_pixel_overlaps,aggregate,NoOverlapError)
+from xagg.auxfuncs import SomeNanWarning
 from xagg.options import set_options
 
 ##### process_weights() tests #####
@@ -607,12 +608,15 @@ def test_aggregate_with_weights(ds=ds):
 	# Get aggregate
 	agg_for = aggregate(ds,wm_for,impl='for_loop')
 	agg_dot = aggregate(ds,wm_dot,impl='dot_product')
+	if _has_numba:
+		agg_num = aggregate(ds,wm_for,impl='numba')
 
 	# Since the "test" for the input ds has [1,8] for the two 
 	# equatorial pixels, the average should just be 4 for timestep 1
 	pd.testing.assert_series_equal(agg_for.agg.test,pd.Series([[[[4,5,6]]]],name='test'))
 	pd.testing.assert_series_equal(agg_dot.agg.test,pd.Series([[[[4,5,6]]]],name='test'))
-
+	if _has_numba:
+		pd.testing.assert_series_equal(agg_num.agg.test,pd.Series([[[[4,5,6]]]],name='test'))
 
 
 def test_aggregate_with_mismatched_grid():
@@ -644,9 +648,13 @@ def test_aggregate_with_mismatched_grid():
 	# Get aggregate
 	agg_for = aggregate(ds,wm_for,impl='for_loop')
 	agg_dot = aggregate(ds,wm_dot,impl='dot_product')
+	if _has_numba:
+		agg_num = aggregate(ds,wm_for,impl='numba')
 
 	pd.testing.assert_series_equal(agg_for.agg.test,pd.Series([[[[18.9999,19.9999,20.9999]]]],name='test'),atol=1e-4)
 	pd.testing.assert_series_equal(agg_dot.agg.test,pd.Series([[[[18.9999,19.9999,20.9999]]]],name='test'),atol=1e-4)
+	if _has_numba:
+		pd.testing.assert_series_equal(agg_num.agg.test,pd.Series([[[[18.9999,19.9999,20.9999]]]],name='test'),atol=1e-4)
 
 
 
@@ -871,6 +879,7 @@ def test_aggregate_with_partialnans():
 						'lon':(['lon'],np.array([0,1])),
 						'bnds':(['bnds'],np.array([0,1])),
 						'time':(['time'],pd.date_range('2019-01-01','2019-01-03'))})
+	# Add some nans in some locations at not all time steps
 	ds['test'] = ds['test'].where(((ds.lon!=1) | (ds.lat!=1) | (ds.time!=pd.to_datetime('2019-01-01'))))
 
 	# get aggregation mapping
@@ -887,12 +896,23 @@ def test_aggregate_with_partialnans():
 	wm_dot = get_pixel_overlaps(gdf,pix_agg,impl='dot_product')
 
 	# Get aggregate
-	with pytest.warns(UserWarning):
+	with pytest.warns(SomeNanWarning):
 		agg_for = aggregate(ds,wm_for,impl='for_loop')
-	with pytest.warns(UserWarning):
+	with pytest.warns(SomeNanWarning):
 		agg_dot = aggregate(ds,wm_dot,impl='dot_product')
 	if _has_numba:
-		agg_num = aggregate(ds,wm_for,impl='numba')
+		with pytest.warns(SomeNanWarning):
+			agg_num = aggregate(ds,wm_for,impl='numba')
+
+#def test_aggregate_with_fullnans():
+	# This is a test to make sure that pixels that have no data / are np.nan
+	# are ignored in the final aggregation, and the remaining overlapping
+	# pixels are renormalized
+
+
+#def test_aggregate_with_polynans():
+	# This is a test to make sure that regions that are completely outside
+	# of the dataset are nan of the same shape as the others have data
 
 ##### aggregate() silencing tests #####
 # Create polygon covering multiple pixels