Clipping large dataset results in out of memory error

[jackieleng]

I'm trying to use rioxarray to clip a large dataset (GeoTIFF file of 5.5GB). However this results in an out of memory error. The GeoTIFF is tiled and I'm using chunks. The clipping does seem to work well for smaller datasets, but is there a way to make it work for larger datasets?

Memory error:

In [1]: %load sess.py

In [2]: # %load sess.py
   ...: import rioxarray
   ...: import xarray as xr
   ...: d = xr.open_rasterio('./data/DSM_Neretva_basin1x1/DSM_Neretva_basin_1x1.tif', chunks=128)
   ...: g = {"type": "Polygon", "coordinates": [[[6485307.1, 4798895.6], [6485728.7, 4799138.2], [6485307.1, 4799138.2], [6485307.1, 4798895.6]]]}
   ...: d.rio.clip(geometries=[g])
   ...:
---------------------------------------------------------------------------
CPLE_OutOfMemoryError                     Traceback (most recent call last)
<ipython-input-2-5021c68037da> in <module>
      4 d = xr.open_rasterio('./data/DSM_Neretva_basin1x1/DSM_Neretva_basin_1x1.tif', chunks=128)
      5 g = {"type": "Polygon", "coordinates": [[[6485307.1, 4798895.6], [6485728.7, 4799138.2], [6485307.1, 4799138.2], [6485307.1, 4798895.6]]]}
----> 6 d.rio.clip(geometries=[g])

~/#/.venv/lib/python3.8/site-packages/rioxarray/rioxarray.py in clip(self, geometries, crs, all_touched, drop, invert)
   1510             ]
   1511
-> 1512         clip_mask_arr = geometry_mask(
   1513             geometries=geometries,
   1514             out_shape=(int(self.height), int(self.width)),

~/#/.venv/lib/python3.8/site-packages/rasterio/env.py in wrapper(*args, **kwds)
    383         else:
    384             with Env.from_defaults():
--> 385                 return f(*args, **kwds)
    386     return wrapper
    387

~/#/.venv/lib/python3.8/site-packages/rasterio/features.py in geometry_mask(geometries, out_shape, transform, all_touched, invert)
     66     fill, mask_value = (0, 1) if invert else (1, 0)
     67
---> 68     return rasterize(
     69         geometries,
     70         out_shape=out_shape,

~/#/.venv/lib/python3.8/site-packages/rasterio/env.py in wrapper(*args, **kwds)
    380     def wrapper(*args, **kwds):
    381         if local._env:
--> 382             return f(*args, **kwds)
    383         else:
    384             with Env.from_defaults():

~/#/.venv/lib/python3.8/site-packages/rasterio/features.py in rasterize(shapes, out_shape, fill, out, transform, all_touched, merge_alg, default_value, dtype)
    355
    356     transform = guard_transform(transform)
--> 357     _rasterize(valid_shapes, out, transform, all_touched, merge_alg)
    358     return out
    359

rasterio/_features.pyx in rasterio._features._rasterize()

rasterio/_io.pyx in rasterio._io.InMemoryRaster.__cinit__()

rasterio/_err.pyx in rasterio._err.exc_wrap_pointer()

CPLE_OutOfMemoryError: memdataset.cpp, 1541: cannot allocate 1x16602363284 bytes

Dataset:

In [5]: d
Out[5]:
<xarray.DataArray (band: 1, y: 122356, x: 135689)>
dask.array<open_rasterio-fe8c0d71fa0f2e35b730362ce6e2bf0f<this-array>, shape=(1, 122356, 135689), dtype=float32, chunksize=(1, 128, 128), chunktype=numpy.ndarray>
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 4.867e+06 4.867e+06 4.867e+06 ... 4.744e+06 4.744e+06
  * x        (x) float64 6.411e+06 6.411e+06 6.411e+06 ... 6.547e+06 6.547e+06
Attributes:
    transform:      (1.0, 0.0, 6411448.162047689, 0.0, -1.0, 4866854.457952311)
    crs:            +proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_...
    res:            (1.0, 1.0)
    is_tiled:       1
    nodatavals:     (-3.4028234663852886e+38,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area

[snowman2]

This will be much simpler when #115 is implemented.

For now, here are alternative options to look into:

1. Clipping with rasterio: Add mask & bbox arguments to open_rasterio #115 (comment)
2. parse_coordinates=False: Add mask & bbox arguments to open_rasterio #115 (comment)

[jackieleng]

Thanks for you help. That is my current workaround to first use clip_box using the polygon bounding box, and afterwards use clip.

Actually I found another way to do this using pure rasterio:

In [7]: import rasterio.mask

In [8]: d = rasterio.open('./data/DSM_Neretva_basin1x1/DSM_Neretva_basin_1x1.tif')

In [9]: rasterio.mask.mask(d, [g], crop=True)
Out[9]:
(array([[[-3.4028235e+38,  5.9820000e+01,  5.1370892e+01, ...,
           9.9109787e+01,  9.9209747e+01, -3.4028235e+38],
         [-3.4028235e+38,  5.9320000e+01,  5.1531410e+01, ...,
          -3.4028235e+38, -3.4028235e+38, -3.4028235e+38],
         [-3.4028235e+38,  5.8910000e+01,  5.1663692e+01, ...,
          -3.4028235e+38, -3.4028235e+38, -3.4028235e+38],
         ...,
         [-3.4028235e+38,  6.1225384e+01,  6.1267384e+01, ...,
          -3.4028235e+38, -3.4028235e+38, -3.4028235e+38],
         [-3.4028235e+38,  6.1217651e+01,  6.1253475e+01, ...,
          -3.4028235e+38, -3.4028235e+38, -3.4028235e+38],
         [-3.4028235e+38,  6.1203785e+01, -3.4028235e+38, ...,
          -3.4028235e+38, -3.4028235e+38, -3.4028235e+38]]], dtype=float32),
 Affine(1.0, 0.0, 6485306.162047689,
        0.0, -1.0, 4799138.457952311))

However, for me this isn't ideal since I would prefer to work with xarray objects only.

[lwasser]

xds.rio.clip_box(minx=xmin, miny=ymin, maxx=xmax, maxy=ymax)

@snowman2 gosh was searching for an alternative to reproject_match to match raster boundaries. Clip_box worked like a charm and barely consumed any resources!! reproject_match is resource intensive and was crashing my hub kernel. I think i understand why. Thank you for this suggestion and for your work with rioxarray as well!!

I have a related question. it seemed to me that chunking the data didn't help much with memory usage for clipping. Is there a good place to better understand how to optimize these types of operations? i have cores i can use but this seemed to be eating up memory.

[snowman2]

I have a related question. it seemed to me that chunking the data didn't help much with memory usage for clipping. Is there a good place to better understand how to optimize these types of operations? i have cores i can use but this seemed to be eating up memory.

This is something that I thing would be worthy of its own topic. Mind opening up another discussion?

[lwasser]

absolutely. after i hit reply i realized i should have done that to begin with :)

[snowman2]

#226