issue_comments: 211045128

In a past life I made side library that wraps rasterio's API to take and return xarray.DataArrays. It provides IO/clip/warp/rasterize operations on DataArrays, which themselves are annotated with the CRS and affine transforms as attributes.

My most common use case was reading disparate rasters and aligning them to the same grid: 1. Use rasterio to load separate spatial rasters over roughly the same area; let's say one is 30-meter satellite and one is 3-meter agricultural yield. Often I'll immediately wrap them in an xarray.DataArray and persist the CRS and affine transform as attributes. 2. Clip the fine-resolution yield array to my area of interest. I can either use rasterio.read(window=(...)) when reading or xarray.sel(x=slice(...), y=slice(...)) post-hoc. 3. I want to overlay the two arrays in a single xarray.Dataset. They may or may not have the same projection but definitely do not the same grid. I'll leave the fine-resolution data untouched to avoid needless resampling, but instead upsample and re-align the coarse-resolution satellite array to match the affine transform of the clipped fine-resolution array. Behind-the-scenes this looks like:

# agriculture and satellite are both 2D DataArrays. satellite_10m = np.zeros(agriculture.shape, dtype=satellite.dtype) rasterio.warp.transform( satellite.values, satellite_10m, src_transform=satellite.attrs['transform'], src_crs=satellite.attrs['crs'], dst_transform=agriculture.attrs['transform'], dst_crs=agriculture.attrs['crs']) # This forces realization of the dataset in memory. I don't do much out-of-core. 4. Stack them into the same Dataset, sharing the transform and coordinate variables from the fine-resolution array. 5. Now use xarray to do cool compuations on the aligned datasets.

Reprojecting or clipping after reading xarray, like I do, goes against @perrygeo's recommendation. So maybe my example is moot, but I really like being able to do this programmatically in python, not CLI.

Even if xarray's new rasterio backend only provides a reader (and not rasterio.warp or rasterio.features functions), Step 3 shows it's very useful if, in addition to the data, the reader will expose the CRS and affine transform objects to the client.

However, if you both expose the transform and realize the coordinate variables, it's possible for them to diverge as the single source of truth. In my above workflow, anytime I clip (step 2) or warp (step 3) data, my side library needed to manually re-set that DataArray's transform and coordinate variables. (This is surely out of scope for rasterio or xarray!)