using reproject_match to resample small GSD to large GSD in small AOI

[joemcglinchy]

Hello!

I have a use case where i'm trying to use reproject_match to resample a dataset with a small GSD (1 meter) to a larger GSD (landsat 30m, sentinel 10m), within a small AOI. I'm leverage the clip operator on both datasets, loaded in full as xarray DataArrays, before calling reproject_match on the small GSD dataset. I've included a repro case, with all files, on google drive here: https://drive.google.com/drive/folders/1O7xAL0uGk799WqZFqsFqGeAJjCnqSYTy?usp=sharing

The Landsat DataArray is quite small, 2x2 pixels. However, the result of reprojecting the smaller GSD dataset to match the Landsat dataset has all values equal to -9999.0, so I'm wondering if there is some kind of mismatch, either in the band coordinate labels or the spatial reference information that is causing trouble behind the scenes. The smaller GSD dataset is coming from an H5 file, so i'm constructing the DataArray 'from scratch', so to speak, and defining its CRS with the write_crs method.

A different problem is surfacing when trying to do the same with the Sentinel2 DataArray. I'm receiving a memory error, which I might expect, but the error traceback says it is trying to allocate an array with dimensions that seem odd to me at first glance. Traceback below.

Can anyone offer thoughts to what may be happening here?

---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-49-8b02410b4012> in <module>
----> 1 temp2 = neon_ds_s2.rio.reproject_match(s2_aoi_ds, Resampling.average)

C:\software\Anaconda3\envs\x-python\lib\site-packages\rioxarray\raster_array.py in reproject_match(self, match_data_array, resampling)
    436             transform=match_data_array.rio.transform(recalc=True),
    437             shape=match_data_array.rio.shape,
--> 438             resampling=resampling,
    439         )
    440 

C:\software\Anaconda3\envs\x-python\lib\site-packages\rioxarray\raster_array.py in reproject(self, dst_crs, resolution, shape, transform, resampling)
    371         )
    372         rasterio.warp.reproject(
--> 373             source=self._obj.values,
    374             destination=dst_data,
    375             src_transform=src_affine,

C:\software\Anaconda3\envs\x-python\lib\site-packages\xarray\core\dataarray.py in values(self)
    567     def values(self) -> np.ndarray:
    568         """The array's data as a numpy.ndarray"""
--> 569         return self.variable.values
    570 
    571     @values.setter

C:\software\Anaconda3\envs\x-python\lib\site-packages\xarray\core\variable.py in values(self)
    508     def values(self):
    509         """The variable's data as a numpy.ndarray"""
--> 510         return _as_array_or_item(self._data)
    511 
    512     @values.setter

C:\software\Anaconda3\envs\x-python\lib\site-packages\xarray\core\variable.py in _as_array_or_item(data)
    270         data = data.get()
    271     else:
--> 272         data = np.asarray(data)
    273     if data.ndim == 0:
    274         if data.dtype.kind == "M":

C:\software\Anaconda3\envs\x-python\lib\site-packages\numpy\core\_asarray.py in asarray(a, dtype, order)
     81 
     82     """
---> 83     return array(a, dtype, copy=False, order=order)
     84 
     85 

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\array\core.py in __array__(self, dtype, **kwargs)
   1374 
   1375     def __array__(self, dtype=None, **kwargs):
-> 1376         x = self.compute()
   1377         if dtype and x.dtype != dtype:
   1378             x = x.astype(dtype)

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\base.py in compute(self, **kwargs)
    165         dask.base.compute
    166         """
--> 167         (result,) = compute(self, traverse=False, **kwargs)
    168         return result
    169 

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\base.py in compute(*args, **kwargs)
    450         postcomputes.append(x.__dask_postcompute__())
    451 
--> 452     results = schedule(dsk, keys, **kwargs)
    453     return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])
    454 

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\threaded.py in get(dsk, result, cache, num_workers, pool, **kwargs)
     82         get_id=_thread_get_id,
     83         pack_exception=pack_exception,
---> 84         **kwargs
     85     )
     86 

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\local.py in get_async(apply_async, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, **kwargs)
    484                         _execute_task(task, data)  # Re-execute locally
    485                     else:
--> 486                         raise_exception(exc, tb)
    487                 res, worker_id = loads(res_info)
    488                 state["cache"][key] = res

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\local.py in reraise(exc, tb)
    314     if exc.__traceback__ is not tb:
    315         raise exc.with_traceback(tb)
--> 316     raise exc
    317 
    318 

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\local.py in execute_task(key, task_info, dumps, loads, get_id, pack_exception)
    220     try:
    221         task, data = loads(task_info)
--> 222         result = _execute_task(task, data)
    223         id = get_id()
    224         result = dumps((result, id))

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\core.py in _execute_task(arg, cache, dsk)
    119         # temporaries by their reference count and can execute certain
    120         # operations in-place.
--> 121         return func(*(_execute_task(a, cache) for a in args))
    122     elif not ishashable(arg):
    123         return arg

C:\software\Anaconda3\envs\x-python\lib\site-packages\dask\array\core.py in getter(a, b, asarray, lock)
     98         lock.acquire()
     99     try:
--> 100         c = a[b]
    101         if asarray:
    102             c = np.asarray(c)

h5py\_objects.pyx in h5py._objects.with_phil.wrapper()

h5py\_objects.pyx in h5py._objects.with_phil.wrapper()

C:\software\Anaconda3\envs\x-python\lib\site-packages\h5py\_hl\dataset.py in __getitem__(self, args)
    560         single_element = selection.mshape == ()
    561         mshape = (1,) if single_element else selection.mshape
--> 562         arr = numpy.ndarray(mshape, new_dtype, order='C')
    563 
    564         # HDF5 has a bug where if the memory shape has a different rank

MemoryError: Unable to allocate 32.0 MiB for an array with shape (256, 256, 256) and data type int16

[snowman2]

I see dask in your errors. Unfortunately reproject and clip do not currently support dask (#119).

[joemcglinchy]

Thanks for the reply, Alan! Were you able to look at the ipython notebook in the google drive folder? I don't see any errors when clipping the data. Is there a way to explicitly tell the reproject_match function to not use dask? I don't have anything in my code explicitly telling any of the modules to use dask.

In the meantime, I'll give that warped VRT workflow a go again, as that has worked for me in the past when I had larger data. However, with the data being so small here, I hadn't expected a similar issue to arise. The gist provided in that issue doesn't account for warping to a new pixel size, so I'll see if I can add that in.