I'll set up a new issue. @Boorhin, I couldn't confirm the weirdness with the small example, but will put in a note to your comment. If you can reproduce the weirdness on the minimal example, would you make a comment to the new issue?

I've made a simpler example of the _FillValue - append issue: python import numpy as np import xarray as xr ds = xr.Dataset({"a": ("x", [3.], {"_FillValue": np.nan})}) m = {} ds.to_zarr(m) ds.to_zarr(m, append_dim="x") raises ValueError: failed to prevent overwriting existing key _FillValue in attrs. This is probably an encoding field used by xarray to describe how a variable is serialized. To proceed, remove this key from the variable's attributes manually. I'd expect this to just work (effectively concatenating the dataset to itself).

The workaround: python m = {} ds.to_zarr(m) del ds.a.attrs["_FillValue"] ds.to_zarr(m, append_dim="x") does the trick, but doesn't look right.

@dcherian, @Boorhin should we make a new (CF-related) issue out of this and try to keep focussing on append and region use-cases here, which seemed to be the initial problem in this thread (probably by going further through your example @Boorhin?).

btw, as a work-around it works when removing the _FillValue from dst.air (you'll likely only want to do this for the append-writes, not the initial write):

python del dst.air.attrs["_FillValue"] dst.to_zarr(m, append_dim="time") works.

But still, this might call for another issue to solve.

Ok, I believe, I've now reproduced your error:

```python import xarray as xr from rasterio.enums import Resampling import numpy as np ds = xr.tutorial.open_dataset('air_temperature').isel(time=0) ds = ds.rio.write_crs('EPSG:4326') dst = ds.rio.reproject('EPSG:3857', shape=(250, 250), resampling=Resampling.bilinear, nodata=np.nan) dst.air.encoding = {} dst = dst.assign(air=dst.air.expand_dims("time"), time=dst.time.expand_dims("time"))

m = {} dst.to_zarr(m) dst.to_zarr(m, append_dim="time") raises: ValueError: failed to prevent overwriting existing key _FillValue in attrs. This is probably an encoding field used by xarray to describe how a variable is serialized. To proceed, remove this key from the variable's attributes manually. ```

This seems to be due to handling of CF-Conventions which might go wrong in the append case: the CFMaskCoder verifies that there isn't any fill value present in the dataset before defining one here. I'd guess in the append case, one wouldn't want to check if the fill value is already defined, but instead want to check that it is the same. However, I don't know a lot about the CF encoding pieces of xarray...

This error ist unrelated to region or append writes. The dataset dst got the _FillValue attribute from rio.reproject ```

dst.air.attrs {... '_FillValue': nan} ```

but still carries encoding-information from ds, i.e.: ```

dst.air.encoding {'source': '...air_temperature.nc', 'original_shape': (2920, 25, 53), 'dtype': dtype('int16'), 'scale_factor': 0.01, 'grid_mapping': 'spatial_ref'} `` The encoding get's picked up byto_zarr, but asnan(the_FillValuefromrio.reproject) can't be expressed as anint16, it's not possible to write that data. You'll have to get rid of the encoding or specify some encoding and_FillValue` which fit together. #5219 might be related.

🤷 can't help any further without a minimal reproducible example here...

I don't yet know a proper answer, but there'd be three observations I have: * The ValueError seems to be related to the handling of CF-Conventions. I don't yet know if that's independent of this issue or if the error only appears in conjunction with this issue. * As far as I understand, appending should be possible without dropping anything (while potentially overwriting some things). * It shouldn't be possible to change _FillValue during appends, because that might require rewriting everything previously written, which you likely don't want. So if _FillValue is different on the append-call, I'd want xarray to produce an error.

While testing a bit further, I found another case which might potentially be dangerous:

```python

ds is the same as above, but chunksize is {"time": 1, "x": 1}

once on the coordinator

ds.to_zarr("test.zarr", compute=False, encoding={"time": {"chunks": [1]}, "x": {"chunks": [1]}})

in parallel

ds.isel(time=slice(0,1), x=slice(0,1)).to_zarr("test.zarr", mode="r+", region={"time": slice(0,1), "x": slice(0,1)}) ds.isel(time=slice(0,1), x=slice(1,2)).to_zarr("test.zarr", mode="r+", region={"time": slice(0,1), "x": slice(1,2)}) ds.isel(time=slice(0,1), x=slice(2,3)).to_zarr("test.zarr", mode="r+", region={"time": slice(0,1), "x": slice(2,3)}) ds.isel(time=slice(1,2), x=slice(0,1)).to_zarr("test.zarr", mode="r+", region={"time": slice(1,2), "x": slice(0,1)}) ds.isel(time=slice(1,2), x=slice(1,2)).to_zarr("test.zarr", mode="r+", region={"time": slice(1,2), "x": slice(1,2)}) ds.isel(time=slice(1,2), x=slice(2,3)).to_zarr("test.zarr", mode="r+", region={"time": slice(1,2), "x": slice(2,3)}) ```

This example doesn't produce any error, but the time and x coordinates are re-written multiple times without any warning. However, I don't yet know how a proper error / warning should be generated in this case. Maybe the check must be if every written variable touches all region-ed dimensions? But maybe thats overly restrictive?

I'm trying to picture some usage scenarios based on incrementally adding timesteps to data on store. I hope these might help to answer questions from above. In particular, I think that append and region options of to_zarr will imply different usage patterns, so might lead to different answers, and mixing terms might lead to confusion.

I'll use the following dataset for demonstration code:

python ds = xr.Dataset({ "T": (("time", "x"), [[1.,2.,3.],[11.,12.,13.]]), }, coords={ "time": (("time",), [21., 22.]), "x": (("x",), [100., 200., 300.]) }).chunk({"time": 1})

append

The purpose of append is to add (one or many) elements along one dimension after the end of all currently existing elements. This implies a read-modify-write cycle to at least the total shape of the array. Furthermore, the place to write new chunks is determined by the current shape of the existing array. Due to these implications, it doesn't seem to be useful to try append in parallel (it would become ambiguous where to write) and it doesn't seem to be too useful (but possible) to only write some of the variables defined on the append-dimension, because all other variables would implicitly be filled with fill_value and those places couldn't be filled with another append anymore.

As a consquence, append-mode writes will always have to be sequential and writes to objects shared touched by multiple append calls will always have a defined behaviour, even if they are modified / overwritten with each call. Creating and appending works as follows:

```python

writes 0-sized time-dimension, so only metadata and non-time dependent variables

ds.isel(time=slice(0,0)).to_zarr("test_append.zarr")

!tree -a test_append.zarr

ds.isel(time=slice(0,1)).to_zarr("test_append.zarr", mode="a", append_dim="time") ds.isel(time=slice(1,2)).to_zarr("test_append.zarr", mode="a", append_dim="time")

print() print("final dataset:") !tree -a test_append.zarr ```

In this case, x would be overwritten with each append call, but the behaviour is well defined as we will only ever append sequentially, so whatever the last write writes into x will be the final result, e.g. [1, 2, 3] in the following case:

python ds.isel(time=slice(0,1)).to_zarr("test_append.zarr", mode="a", append_dim="time") ds2 = ds.assign(x=[1,2,3]) ds2.isel(time=slice(1,2)).to_zarr("test_append.zarr", mode="a", append_dim="time")

If instead, x shouldn't be overwritten, it's possible to append using: python ds.drop(["x"]).isel(time=slice(0,1)).to_zarr("test_append.zarr", mode="a", append_dim="time") ds.drop(["x"]).isel(time=slice(1,2)).to_zarr("test_append.zarr", mode="a", append_dim="time") This also works already with current xarray and has well defined behaviour. However, if there are many time-independent variables, it might be easier if something like .drop_if_not("time") or something similar would be available.

region

region behaves quite differently from append. It does not modify the shape of the arrays and it does not depend on the shape's value to determine where to write new data (it requires user input to do so). This generally enables parallel writes to the same dataset (if only distinct chunks are touched). But as metadata (e.g. shape) is still shared, updates to metadata must be done in a coordinated (likely sequential) manner.

Generally, the workflow with region would imply writing the metadata once and maybe update it from time to time but sequentially (e.g. on a coordinating node) and write all the chunks in parallel on worker nodes, while carefully ensuring that no common chunks are overwritten. Let's see how this might look like:

```python ds.to_zarr("test.zarr", compute=False, encoding={"time": {"chunks": [1]}}) !rm test.zarr/time/0 !rm test.zarr/time/1

!tree -a test.zarr

NOTE: these may run in parallel (even if that's not useful in time, but region might also be in time and space)

ds.drop(['x']).isel(time=slice(0,1)).to_zarr("test.zarr", mode="r+", region={"time": slice(0,1)}) ds.drop(['x']).isel(time=slice(1,2)).to_zarr("test.zarr", mode="r+", region={"time": slice(1,2)})

print() print("final dataset:") !tree -a test.zarr ```

The above works and as far as I understand does what we'd want for parallel writes. It also avoids the mentioned ambiguous cases (due to the drop(['x']) statements). However this case is even more cumbersome to write than in the append case. The parallel writes might benefit from again from something like .drop_if_not("time") (which probably can't be optional in this case due to ambiguity). But what's even more problematic is the initial write of array metadata. In order to start building the dataset, I'll have to scaffold an (potentially not yet computed) Dataset of full size and use compute=False to write only metadata. However, this fails for coordinate variables (like time), because those are eagerly loaded and will still be written out. That's why I've removed those chunks in the example above.

If region should be used for parallel append, then there must be some process on a coordinating node which updates the metadata keys (at least by increasing the shape). I don't yet see how that could be written nicely using xarray.

So based on these two kinds of tasks, it seems to me that the actual append and region write-modes of to_zarr are already doing what they should do, but there could be some more convenience functions which would make those tasks much simpler:

• some method like drop_if_not (maybe with a better name) which drops all the things we don't want to keep (maybe we should call it keep instead of drop). This method would essentially result in and simplify mode 1 in @shoyer's answer, which I'd argue is what we actually want in both use cases, becasue the dropped data would already have been written by the coordinating process. I'd believe that mode 1 shouldn't be the default for to_zarr though, because silently dropping data from being written isn't nice to the user.
• some better tooling for writing and updating zarr dataset metadata (I don't know if that would fit in the realm of xarray though, as it looks like handling Datasets without content. For "appending" metadata, I really don't know how I'd picture this propery in xarray world.)

You are right, the coordinates should not be dropped.

I think the function _validate_region has a bug. Currently it checks for all ds.variables if at least one of their dimensions agrees with the ones given in the region argument. However, ds.variables also returns the coordinates. However, we actually only want to check if the ds.data_vars have a dimension intersecting with the given region.

Changing the function to `python def _validate_region(ds, region): if not isinstance(region, dict): raise TypeError(f"region`` must be a dict, got {type(region)}")

for k, v in region.items():
    if k not in ds.dims:
        raise ValueError(
            f"all keys in ``region`` are not in Dataset dimensions, got "
            f"{list(region)} and {list(ds.dims)}"
        )
    if not isinstance(v, slice):
        raise TypeError(
            "all values in ``region`` must be slice objects, got "
            f"region={region}"
        )
    if v.step not in {1, None}:
        raise ValueError(
            "step on all slices in ``region`` must be 1 or None, got "
            f"region={region}"
        )

non_matching_vars = [
    k for k, v in ds.data_vars.items() if not set(region).intersection(v.dims)
]
if non_matching_vars:
    raise ValueError(
        f"when setting `region` explicitly in to_zarr(), all "
        f"variables in the dataset to write must have at least "
        f"one dimension in common with the region's dimensions "
        f"{list(region.keys())}, but that is not "
        f"the case for some variables here. To drop these variables "
        f"from this dataset before exporting to zarr, write: "
        f".drop({non_matching_vars!r})"
    )

``` seems to work.

Hi @Boorhin, I just ran into the same issue. The region argument has to be of type slice, in your case slice(t) instead of just t works:

python import xarray as xr from datetime import datetime,timedelta import numpy as np dt= datetime.now() times= np.arange(dt,dt+timedelta(days=6), timedelta(hours=1)) nodesx,nodesy,layers=np.arange(10,50), np.arange(10,50)+15, np.arange(10) ds=xr.Dataset() ds.coords['time']=('time', times) ds.coords['node_x']=('node', nodesx) ds.coords['node_y']=('node', nodesy) ds.coords['layer']=('layer', layers) outfile='my_zarr' varnames=['potato','banana', 'apple'] for var in varnames: ds[var]=(('time', 'layer', 'node'), np.zeros((len(times), len(layers),len(nodesx)))) ds.to_zarr(outfile, mode='a') for t in range(len(times)): for var in varnames: ds[var].isel(time=slice(t)).values += np.random.random((len(layers),len(nodesx))) ds.isel(time=slice(t)).to_zarr(outfile, region={"time": slice(t)})

This leads however to another issue: ```python

ValueError Traceback (most recent call last) <ipython-input-52-bb3d2c1adc12> in <module> 18 for var in varnames: 19 ds[var].isel(time=slice(t)).values += np.random.random((len(layers),len(nodesx))) ---> 20 ds.isel(time=slice(t)).to_zarr(outfile, region={"time": slice(t)})

~/.local/lib/python3.8/site-packages/xarray/core/dataset.py in to_zarr(self, store, chunk_store, mode, synchronizer, group, encoding, compute, consolidated, append_dim, region, safe_chunks) 2029 encoding = {} 2030 -> 2031 return to_zarr( 2032 self, 2033 store=store,

~/.local/lib/python3.8/site-packages/xarray/backends/api.py in to_zarr(dataset, store, chunk_store, mode, synchronizer, group, encoding, compute, consolidated, append_dim, region, safe_chunks) 1359 1360 if region is not None: -> 1361 _validate_region(dataset, region) 1362 if append_dim is not None and append_dim in region: 1363 raise ValueError(

~/.local/lib/python3.8/site-packages/xarray/backends/api.py in _validate_region(ds, region) 1272 ] 1273 if non_matching_vars: -> 1274 raise ValueError( 1275 f"when setting region explicitly in to_zarr(), all " 1276 f"variables in the dataset to write must have at least "

ValueError: when setting region explicitly in to_zarr(), all variables in the dataset to write must have at least one dimension in common with the region's dimensions ['time'], but that is not the case for some variables here. To drop these variables from this dataset before exporting to zarr, write: .drop(['node_x', 'node_y', 'layer']) ```

Here, the solution is however provided with the error message. Following the instructions, the snippet below finally works (as far as I can tell):

```python import xarray as xr from datetime import datetime,timedelta import numpy as np dt= datetime.now() times= np.arange(dt,dt+timedelta(days=6), timedelta(hours=1)) nodesx,nodesy,layers=np.arange(10,50), np.arange(10,50)+15, np.arange(10) ds=xr.Dataset() ds.coords['time']=('time', times)

ds.coords['node_x']=('node', nodesx)

ds.coords['node_y']=('node', nodesy)

ds.coords['layer']=('layer', layers)

outfile='my_zarr' varnames=['potato','banana', 'apple'] for var in varnames: ds[var]=(('time', 'layer', 'node'), np.zeros((len(times), len(layers),len(nodesx)))) ds.to_zarr(outfile, mode='a') for t in range(len(times)): for var in varnames: ds[var].isel(time=slice(t)).values += np.random.random((len(layers),len(nodesx))) ds.isel(time=slice(t)).to_zarr(outfile, region={"time": slice(t)}) ```

Maybe one would like to generalise region in api.py to allow for single indices or throw a hint in case an a type different to a slice is provided.

Cheers