Add option/ability to use Rasterio MemoryFiles as outputs/destinations

[ciaransweet]

Hey folks,

I've been a heavy rasterio user for a while, working on a project currently that is going to be taking NetCDF4 files and converting them to COGs (with some other stuff in between).

I found rioxarray really handy for getting my NetCDF4 files into Rasters, what would be really awesome is integration/the ability to use rasterio MemoryFiles.

The use case being that really I'd rather not have to write to disk and clean up intermediate files (I do a lot of serverless based workloads).

For example I could see my code looking something like:

from rasterio.io import MemoryFile
import rioxarray
import xarray

xds = xarray.open("my-dataset.nc4")
# Various rioxarray stuff
with MemoryFile() as memfile:
    with memfile.open(driver='GTiff', count=3, ...) as dataset:
        dataset.write(xds.rio.to_dataset()) # New function?
        # or
        dataset.write(xds.rio.to_raster()) # By default a Rasterio Dataset is returned if no filepath provided?

This would then allow me to work on the resultant raster within Rasterio without worrying about writing to disk and cleaning up intermediate files

Hopefully I've explained myself well, happy to try and flesh it out if needed :D

Cheers!

[Zeitsperre]

Neat idea. I'd love to see this as well!

[snowman2]

This would then allow me to work on the resultant raster within Rasterio without worrying about writing to disk and cleaning up intermediate files

Curious what you need to do in rasterio that you cannot do in rioxarray with xds.copy()?

Way to do it in current state

This seems to work:

from rasterio.io import MemoryFile
import xarray

xds = xarray.open("my-dataset.nc4")

with MemoryFile() as memfile:
    xds.rio.to_raster(memfile.name)
    with memfile.open() as dataset:
        ....

Potential future implementation

I think the logic in to_raster could be re-used and could simplify some things and make it more consistent:

with xds.rio.to_memory_file() as memfile:
    ...

[ciaransweet]

I have other libraries that are plugins for rasterio, hence the requirement to be able to take them to rasterio memory files.

I'll be honest I'd not thought about using .values (I wasn't really aware) - That probably would work yes!

I do like your suggested future implementation though also :D

[ciaransweet]

Hey @snowman2, I just gave that a go and it worked a charm. No rush in the memory file implementation 😄

Many thanks for your help!

[RichardScottOZ]

As far as rasterio goes, not something I have tried - writing out many band rasters.

E.g. I have plenty of 20+ band geotiffs lying around and working on a 27 band one today, the example I think just mention the output of one band?

[RichardScottOZ]

Thanks - I remembered reading this in to_raster :-

**profile_kwargs – Additional keyword arguments to pass into writing the raster. The nodata, transform, crs, count, width, and height attributes are ignored.

So do you have a quick example of such to point to, please?

[snowman2]

An example is in to_raster

rioxarray/rioxarray/raster_array.py

Lines 870 to 917 in 29fb395

	if driver is None and LooseVersion(rasterio.__version__) < LooseVersion("1.2"):
	driver = "GTiff"
	dtype = str(self._obj.dtype) if dtype is None else dtype
	# get the output profile from the rasterio object
	# if opened with xarray.open_rasterio()
	try:
	out_profile = self._obj._file_obj.acquire().profile
	except AttributeError:
	out_profile = {}
	out_profile.update(profile_kwargs)
	# filter out the generated attributes
	out_profile = {
	key: value
	for key, value in out_profile.items()
	if key
	not in (
	"driver",
	"height",
	"width",
	"crs",
	"transform",
	"nodata",
	"count",
	"dtype",
	)
	}
	rio_nodata = (
	self.encoded_nodata if self.encoded_nodata is not None else self.nodata
	)
	if rio_nodata is not None:
	# Ensure dtype of output data matches the expected dtype.
	# This check is added here as the dtype of the data is
	# converted right before writing.
	rio_nodata = _ensure_nodata_dtype(rio_nodata, dtype)
	return RasterioWriter(raster_path=raster_path).to_raster(
	xarray_dataarray=self._obj,
	tags=tags,
	driver=driver,
	height=int(self.height),
	width=int(self.width),
	count=int(self.count),
	dtype=dtype,
	crs=self.crs,
	transform=self.transform(recalc=recalc_transform),
	nodata=rio_nodata,

[RichardScottOZ]

Yes, maybe I am misunderstanding how that is working?- isn't that explicitly avoiding setting count? e.g. if you change it, will be skipped?

That is the code for how it works, not this is 'getting a metadata dictionary from rasterio and updating that' and passing as you would using that type example.

[snowman2]

The attributes on the dataarray are what is used to generate the profile in to_raster. It does not allow you to overwrite them when writing to a raster as the attributes on the dataarray are the only ones that should be used.

So, if you wanted to use the same thing for the memfile:

with MemoryFile() as memfile:
    with memfile.open(
        driver='GTiff',
        count=int(rds.rio.count),
        width=int(rds.rio.width)
        height=int(rds.rio.height),
        dtype=str(rds.dtype),
        crs=rds.rio.crs, 
        transform=rds.rio.transform(recalc=True), 
        nodata=rds.rio.nodata, 
        **rasterio_kwargs
    ) as dataset:

[RichardScottOZ]

Thanks for explaining!

That is really good.

I'll see if I can do an example later that illustrates the attributes.

[RichardScottOZ]

Odd one here, I get

---------------------------------------------------------------------------
InvalidDimensionOrder                     Traceback (most recent call last)
<ipython-input-236-49a27e007153> in <module>
      2 rds = result
      3 with MemoryFile() as memfile:
----> 4     with memfile.open(driver='GTiff',count=int(rds.rio.count),width=int(rds.rio.width),height=int(rds.rio.height),dtype=str(rds.dtype),crs=rds.rio.crs, transform=rds.rio.transform(recalc=True), nodata=rds.rio.nodata 
      5     ) as dataset:
      6         data = rds.values

/env/lib/python3.8/site-packages/rioxarray/rioxarray.py in count(self)
    797         if self._count is not None:
    798             return self._count
--> 799         extra_dim = self._check_dimensions()
    800         self._count = 1
    801         if extra_dim is not None:

/env/lib/python3.8/site-packages/rioxarray/rioxarray.py in _check_dimensions(self)
    784             )
    785         if not extra_dims and self._obj.dims != (self.y_dim, self.x_dim):
--> 786             raise InvalidDimensionOrder(
    787                 "Invalid dimension order. Expected order: {0}"
    788                 "You can use `DataArray.transpose{0}` "

InvalidDimensionOrder: Invalid dimension order. Expected order: ('y', 'x')You can use `DataArray.transpose('y', 'x')` to reorder your dimensions.

[RichardScottOZ]

[RichardScottOZ]

rds = result.to_array()

[snowman2]

@RichardScottOZ if you think there is a bug, please open up an issue with more details.

[RichardScottOZ]

Not sure yet, but making it an array works, anyway. e.g. if I am erroneously assuming can just write a dataset?

Other than that, have been putting your clarification above to good work, so thanks again!