Xarray supports direct serialization and IO to several file formats, fromsimple Pickle files to the more flexible netCDFformat (recommended).
netCDF#
The recommended way to store xarray data structures is netCDF, whichis a binary file format for self-described datasets that originatedin the geosciences. Xarray is based on the netCDF data model, so netCDF fileson disk directly correspond to Dataset objects (more accurately,a group in a netCDF file directly corresponds to a Dataset object.See Groups for more.)
NetCDF is supported on almost all platforms, and parsers existfor the vast majority of scientific programming languages. Recent versions ofnetCDF are based on the even more widely used HDF5 file-format.
If you aren’t familiar with this data format, the netCDF FAQ is a goodplace to start.
Reading and writing netCDF files with xarray requires scipy, h5netcdf, or thenetCDF4-Python library to be installed. SciPy only supports reading and writingof netCDF V3 files.
We can save a Dataset to disk using theDataset.to_netcdf() method:
In [1]: ds = xr.Dataset( ...: {"foo": (("x", "y"), np.random.rand(4, 5))}, ...: coords={ ...: "x": [10, 20, 30, 40], ...: "y": pd.date_range("2000-01-01", periods=5), ...: "z": ("x", list("abcd")), ...: }, ...: ) ...: In [2]: ds.to_netcdf("saved_on_disk.nc")
By default, the file is saved as netCDF4 (assuming netCDF4-Python isinstalled). You can control the format and engine used to write the file withthe format and engine arguments.
Tip
Using the h5netcdf packageby passing engine='h5netcdf' to open_dataset() cansometimes be quicker than the default engine='netcdf4' that uses thenetCDF4 package.
We can load netCDF files to create a new Dataset usingopen_dataset():
In [3]: ds_disk = xr.open_dataset("saved_on_disk.nc")In [4]: ds_diskOut[4]: <xarray.Dataset> Size: 248BDimensions: (x: 4, y: 5)Coordinates: * x (x) int64 32B 10 20 30 40 * y (y) datetime64[ns] 40B 2000-01-01 2000-01-02 ... 2000-01-05 z (x) <U1 16B ...Data variables: foo (x, y) float64 160B ...
Similarly, a DataArray can be saved to disk using theDataArray.to_netcdf() method, and loadedfrom disk using the open_dataarray() function. As netCDF filescorrespond to Dataset objects, these functions internallyconvert the DataArray to a Dataset before saving, and then convert backwhen loading, ensuring that the DataArray that is loaded is always exactlythe same as the one that was saved.
A dataset can also be loaded or written to a specific group within a netCDFfile. To load from a group, pass a group keyword argument to theopen_dataset function. The group can be specified as a path-likestring, e.g., to access subgroup ‘bar’ within group ‘foo’ pass‘/foo/bar’ as the group argument. When writing multiple groups in one file,pass mode='a' to to_netcdf to ensure that each call does not delete thefile.
Data is always loaded lazily from netCDF files. You can manipulate, slice and subsetDataset and DataArray objects, and no array values are loaded into memory untilyou try to perform some sort of actual computation. For an example of how theselazy arrays work, see the OPeNDAP section below.
There may be minor differences in the Dataset object returnedwhen reading a NetCDF file with different engines.
It is important to note that when you modify values of a Dataset, even onelinked to files on disk, only the in-memory copy you are manipulating in xarrayis modified: the original file on disk is never touched.
Tip
Xarray’s lazy loading of remote or on-disk datasets is often but not alwaysdesirable. Before performing computationally intense operations, it isoften a good idea to load a Dataset (or DataArray) entirely into memory byinvoking the Dataset.load() method.
Datasets have a Dataset.close() method to close the associatednetCDF file. However, it’s often cleaner to use a with statement:
# this automatically closes the dataset after useIn [5]: with xr.open_dataset("saved_on_disk.nc") as ds: ...: print(ds.keys()) ...: KeysView(<xarray.Dataset> Size: 248BDimensions: (x: 4, y: 5)Coordinates: * x (x) int64 32B 10 20 30 40 * y (y) datetime64[ns] 40B 2000-01-01 2000-01-02 ... 2000-01-05 z (x) <U1 16B ...Data variables: foo (x, y) float64 160B ...)
Although xarray provides reasonable support for incremental reads of files ondisk, it does not support incremental writes, which can be a useful strategyfor dealing with datasets too big to fit into memory. Instead, xarray integrateswith dask.array (see Parallel computing with Dask), which provides a fully featured engine forstreaming computation.
It is possible to append or overwrite netCDF variables using the mode='a'argument. When using this option, all variables in the dataset will be writtento the original netCDF file, regardless if they exist in the original dataset.
Groups#
NetCDF groups are not supported as part of the Dataset data model.Instead, groups can be loaded individually as Dataset objects.To do so, pass a group keyword argument to theopen_dataset() function. The group can be specified as a path-likestring, e.g., to access subgroup 'bar' within group 'foo' pass'/foo/bar' as the group argument.
In a similar way, the group keyword argument can be given to theDataset.to_netcdf() method to write to a groupin a netCDF file.When writing multiple groups in one file, pass mode='a' toDataset.to_netcdf() to ensure that each call does not delete the file.For example:
In [6]: ds1 = xr.Dataset({"a": 0})In [7]: ds2 = xr.Dataset({"b": 1})In [8]: ds1.to_netcdf("file.nc", group="A")In [9]: ds2.to_netcdf("file.nc", group="B", mode="a")
We can verify that two groups have been saved using the ncdump command-line utility.
$ ncdump file.ncnetcdf file {group: A { variables: int64 a ; data: a = 0 ; } // group Agroup: B { variables: int64 b ; data: b = 1 ; } // group B}
Either of these groups can be loaded from the file as an independent Dataset object:
In [10]: group1 = xr.open_dataset("file.nc", group="A")In [11]: group1Out[11]: <xarray.Dataset>Dimensions: ()Data variables: a int64 ...In [12]: group2 = xr.open_dataset("file.nc", group="B")In [13]: group2Out[13]: <xarray.Dataset>Dimensions: ()Data variables: b int64 ...
Note
For native handling of multiple groups with xarray, including I/O, you might be interested in the experimentalxarray-datatree package.
Reading encoded data#
NetCDF files follow some conventions for encoding datetime arrays (as numberswith a “units” attribute) and for packing and unpacking data (asdescribed by the “scale_factor” and “add_offset” attributes). If the argumentdecode_cf=True (default) is given to open_dataset(), xarray will attemptto automatically decode the values in the netCDF objects according toCF conventions. Sometimes this will fail, for example, if a variablehas an invalid “units” or “calendar” attribute. For these cases, you canturn this decoding off manually.
You can view this encoding information (among others) in theDataArray.encoding andDataArray.encoding attributes:
In [14]: ds_disk["y"].encodingOut[14]: {'dtype': dtype('int64'), 'zlib': False, 'szip': False, 'zstd': False, 'bzip2': False, 'blosc': False, 'shuffle': False, 'complevel': 0, 'fletcher32': False, 'contiguous': True, 'chunksizes': None, 'source': '/home/docs/checkouts/readthedocs.org/user_builds/xray/checkouts/stable/doc/saved_on_disk.nc', 'original_shape': (5,), 'units': 'days since 2000-01-01 00:00:00', 'calendar': 'proleptic_gregorian'}In [15]: ds_disk.encodingOut[15]: {'unlimited_dims': set(), 'source': '/home/docs/checkouts/readthedocs.org/user_builds/xray/checkouts/stable/doc/saved_on_disk.nc'}
Note that all operations that manipulate variables other than indexingwill remove encoding information.
In some cases it is useful to intentionally reset a dataset’s original encoding values.This can be done with either the Dataset.drop_encoding() orDataArray.drop_encoding() methods.
In [16]: ds_no_encoding = ds_disk.drop_encoding()In [17]: ds_no_encoding.encodingOut[17]: {}
Reading multi-file datasets#
NetCDF files are often encountered in collections, e.g., with different filescorresponding to different model runs or one file per timestamp.Xarray can straightforwardly combine such files into a single Dataset by making use ofconcat(), merge(), combine_nested() andcombine_by_coords(). For details on the difference between thesefunctions see Combining data.
Xarray includes support for manipulating datasets that don’t fit into memorywith dask. If you have dask installed, you can open multiple filessimultaneously in parallel using open_mfdataset():
xr.open_mfdataset('my/files/*.nc', parallel=True)
This function automatically concatenates and merges multiple files into asingle xarray dataset.It is the recommended way to open multiple files with xarray.For more details on parallel reading, see Combining along multiple dimensions, Reading and writing data and ablog post by Stephan Hoyer.open_mfdataset() takes many kwargs that allow you tocontrol its behaviour (for e.g. parallel, combine, compat, join, concat_dim).See its docstring for more details.
Note
A common use-case involves a dataset distributed across a large number of files witheach file containing a large number of variables. Commonly, a few of these variablesneed to be concatenated along a dimension (say "time"), while the rest are equalacross the datasets (ignoring floating point differences). The following commandwith suitable modifications (such as parallel=True) works well with such datasets:
xr.open_mfdataset('my/files/*.nc', concat_dim="time", combine="nested", data_vars='minimal', coords='minimal', compat='override')
This command concatenates variables along the "time" dimension, but only those thatalready contain the "time" dimension (data_vars='minimal', coords='minimal').Variables that lack the "time" dimension are taken from the first dataset(compat='override').
Sometimes multi-file datasets are not conveniently organized for easy use of open_mfdataset().One can use the preprocess argument to provide a function that takes a datasetand returns a modified Dataset.open_mfdataset() will call preprocess on every dataset(corresponding to each file) prior to combining them.
If open_mfdataset() does not meet your needs, other approaches are possible.The general pattern for parallel reading of multiple filesusing dask, modifying those datasets and then combining into a single Dataset is:
def modify(ds): # modify ds here return ds# this is basically what open_mfdataset doesopen_kwargs = dict(decode_cf=True, decode_times=False)open_tasks = [dask.delayed(xr.open_dataset)(f, **open_kwargs) for f in file_names]tasks = [dask.delayed(modify)(task) for task in open_tasks]datasets = dask.compute(tasks) # get a list of xarray.Datasetscombined = xr.combine_nested(datasets) # or some combination of concat, merge
As an example, here’s how we could approximate MFDataset from the netCDF4library:
from glob import globimport xarray as xrdef read_netcdfs(files, dim): # glob expands paths with * to a list of files, like the unix shell paths = sorted(glob(files)) datasets = [xr.open_dataset(p) for p in paths] combined = xr.concat(datasets, dim) return combinedcombined = read_netcdfs('/all/my/files/*.nc', dim='time')
This function will work in many cases, but it’s not very robust. First, itnever closes files, which means it will fail if you need to load more thana few thousand files. Second, it assumes that you want all the data from eachfile and that it can all fit into memory. In many situations, you only needa small subset or an aggregated summary of the data from each file.
Here’s a slightly more sophisticated example of how to remedy thesedeficiencies:
def read_netcdfs(files, dim, transform_func=None): def process_one_path(path): # use a context manager, to ensure the file gets closed after use with xr.open_dataset(path) as ds: # transform_func should do some sort of selection or # aggregation if transform_func is not None: ds = transform_func(ds) # load all data from the transformed dataset, to ensure we can # use it after closing each original file ds.load() return ds paths = sorted(glob(files)) datasets = [process_one_path(p) for p in paths] combined = xr.concat(datasets, dim) return combined# here we suppose we only care about the combined mean of each file;# you might also use indexing operations like .sel to subset datasetscombined = read_netcdfs('/all/my/files/*.nc', dim='time', transform_func=lambda ds: ds.mean())
This pattern works well and is very robust. We’ve used similar code to processtens of thousands of files constituting 100s of GB of data.
Writing encoded data#
Conversely, you can customize how xarray writes netCDF files on disk byproviding explicit encodings for each dataset variable. The encodingargument takes a dictionary with variable names as keys and variable specificencodings as values. These encodings are saved as attributes on the netCDFvariables on disk, which allows xarray to faithfully read encoded data back intomemory.
It is important to note that using encodings is entirely optional: if you do notsupply any of these encoding options, xarray will write data to disk using adefault encoding, or the options in the encoding attribute, if set.This works perfectly fine in most cases, but encoding can be useful foradditional control, especially for enabling compression.
In the file on disk, these encodings are saved as attributes on each variable, whichallow xarray and other CF-compliant tools for working with netCDF files to correctlyread the data.
Scaling and type conversions#
These encoding options work on any version of the netCDF file format:
dtype: Any valid NumPy dtype or string convertible to a dtype, e.g.,'int16'or'float32'. This controls the type of the data written on disk._FillValue: Values ofNaNin xarray variables are remapped to this value whensaved on disk. This is important when converting floating point with missing valuesto integers on disk, becauseNaNis not a valid value for integer dtypes. Bydefault, variables with float types are attributed a_FillValueofNaNin theoutput file, unless explicitly disabled with an encoding{'_FillValue': None}.scale_factorandadd_offset: Used to convert from encoded data on disk toto the decoded data in memory, according to the formuladecoded = scale_factor * encoded + add_offset.
These parameters can be fruitfully combined to compress discretized data on disk. Forexample, to save the variable foo with a precision of 0.1 in 16-bit integers whileconverting NaN to -9999, we would useencoding={'foo': {'dtype': 'int16', 'scale_factor': 0.1, '_FillValue': -9999}}.Compression and decompression with such discretization is extremely fast.
String encoding#
Xarray can write unicode strings to netCDF files in two ways:
As variable length strings. This is only supported on netCDF4 (HDF5) files.
By encoding strings into bytes, and writing encoded bytes as a characterarray. The default encoding is UTF-8.
By default, we use variable length strings for compatible files and fall-backto using encoded character arrays. Character arrays can be selected even fornetCDF4 files by setting the dtype field in encoding to S1(corresponding to NumPy’s single-character bytes dtype).
If character arrays are used:
The string encoding that was used is stored ondisk in the
_Encodingattribute, which matches an ad-hoc conventionadopted by the netCDF4-Python library.At the time of this writing (October 2017), a standard convention for indicatingstring encoding for character arrays in netCDF files wasstill under discussion.Technically, you can useany string encoding recognized by Python if you feel the need to deviate from UTF-8,by setting the_Encodingfield inencoding. Butwe don’t recommend it.The character dimension name can be specified by the
char_dim_namefield of a variable’sencoding. If the name of the character dimension is not specified, the default isf'string{data.shape[-1]}'. When decoding character arrays from existing files, thechar_dim_nameis added to the variablesencodingto preserve if encoding happens, butthe field can be edited by the user.
Warning
Missing values in bytes or unicode string arrays (represented by NaN inxarray) are currently written to disk as empty strings ''. This meansmissing values will not be restored when data is loaded from disk.This behavior is likely to change in the future (GH1647).Unfortunately, explicitly setting a _FillValue for string arrays to handlemissing values doesn’t work yet either, though we also hope to fix this in thefuture.
Chunk based compression#
zlib, complevel, fletcher32, contiguous and chunksizescan be used for enabling netCDF4/HDF5’s chunk based compression, as describedin the documentation for createVariable for netCDF4-Python. This only worksfor netCDF4 files and thus requires using format='netCDF4' and eitherengine='netcdf4' or engine='h5netcdf'.
Chunk based gzip compression can yield impressive space savings, especiallyfor sparse data, but it comes with significant performance overhead. HDF5libraries can only read complete chunks back into memory, and maximumdecompression speed is in the range of 50-100 MB/s. Worse, HDF5’s compressionand decompression currently cannot be parallelized with dask. For these reasons, werecommend trying discretization based compression (described above) first.
Time units#
The units and calendar attributes control how xarray serializes datetime64 andtimedelta64 arrays to datasets on disk as numeric values. The units encodingshould be a string like 'days since 1900-01-01' for datetime64 data or a stringlike 'days' for timedelta64 data. calendar should be one of the calendar typessupported by netCDF4-python: ‘standard’, ‘gregorian’, ‘proleptic_gregorian’ ‘noleap’,‘365_day’, ‘360_day’, ‘julian’, ‘all_leap’, ‘366_day’.
By default, xarray uses the 'proleptic_gregorian' calendar and units of the smallest timedifference between values, with a reference time of the first time value.
Coordinates#
You can control the coordinates attribute written to disk by specifying DataArray.encoding["coordinates"].If not specified, xarray automatically sets DataArray.encoding["coordinates"] to a space-delimited listof names of coordinate variables that share dimensions with the DataArray being written.This allows perfect roundtripping of xarray datasets but may not be desirable.When an xarray Dataset contains non-dimensional coordinates that do not share dimensions with any ofthe variables, these coordinate variable names are saved under a “global” "coordinates" attribute.This is not CF-compliant but again facilitates roundtripping of xarray datasets.
Invalid netCDF files#
The library h5netcdf allows writing some dtypes (booleans, complex, …) that aren’tallowed in netCDF4 (seeh5netcdf documentation).This feature is available through DataArray.to_netcdf() andDataset.to_netcdf() when used with engine="h5netcdf"and currently raises a warning unless invalid_netcdf=True is set:
# Writing complex valued dataIn [18]: da = xr.DataArray([1.0 + 1.0j, 2.0 + 2.0j, 3.0 + 3.0j])In [19]: da.to_netcdf("complex.nc", engine="h5netcdf", invalid_netcdf=True)# Reading it backIn [20]: reopened = xr.open_dataarray("complex.nc", engine="h5netcdf")In [21]: reopenedOut[21]: <xarray.DataArray (dim_0: 3)> Size: 48B[3 values with dtype=complex128]Dimensions without coordinates: dim_0
Warning
Note that this produces a file that is likely to be not readable by other netCDFlibraries!
HDF5#
HDF5 is both a file format and a data model for storing information. HDF5 storesdata hierarchically, using groups to create a nested structure. HDF5 is a moregeneral version of the netCDF4 data model, so the nested structure is one of manysimilarities between the two data formats.
Reading HDF5 files in xarray requires the h5netcdf engine, which can be installedwith conda install h5netcdf. Once installed we can use xarray to open HDF5 files:
xr.open_dataset("/path/to/my/file.h5")
The similarities between HDF5 and netCDF4 mean that HDF5 data can be written with thesame Dataset.to_netcdf() method as used for netCDF4 data:
In [22]: ds = xr.Dataset( ....: {"foo": (("x", "y"), np.random.rand(4, 5))}, ....: coords={ ....: "x": [10, 20, 30, 40], ....: "y": pd.date_range("2000-01-01", periods=5), ....: "z": ("x", list("abcd")), ....: }, ....: ) ....: In [23]: ds.to_netcdf("saved_on_disk.h5")
Groups#
If you have multiple or highly nested groups, xarray by default may not read the groupthat you want. A particular group of an HDF5 file can be specified using the groupargument:
xr.open_dataset("/path/to/my/file.h5", group="/my/group")
While xarray cannot interrogate an HDF5 file to determine which groups are available,the HDF5 Python reader h5py can be used instead.
Natively the xarray data structures can only handle one level of nesting, organized asDataArrays inside of Datasets. If your HDF5 file has additional levels of hierarchy youcan only access one group and a time and will need to specify group names.
Note
For native handling of multiple HDF5 groups with xarray, including I/O, you might beinterested in the experimentalxarray-datatree package.
Zarr#
Zarr is a Python package that provides an implementation of chunked, compressed,N-dimensional arrays.Zarr has the ability to store arrays in a range of ways, including in memory,in files, and in cloud-based object storage such as Amazon S3 andGoogle Cloud Storage.Xarray’s Zarr backend allows xarray to leverage these capabilities, includingthe ability to store and analyze datasets far too large fit onto disk(particularly in combination with dask).
Xarray can’t open just any zarr dataset, because xarray requires specialmetadata (attributes) describing the dataset dimensions and coordinates.At this time, xarray can only open zarr datasets with these special attributes,such as zarr datasets written by xarray,netCDF,or GDAL.For implementation details, see Zarr Encoding Specification.
To write a dataset with zarr, we use the Dataset.to_zarr() method.
To write to a local directory, we pass a path to a directory:
In [24]: ds = xr.Dataset( ....: {"foo": (("x", "y"), np.random.rand(4, 5))}, ....: coords={ ....: "x": [10, 20, 30, 40], ....: "y": pd.date_range("2000-01-01", periods=5), ....: "z": ("x", list("abcd")), ....: }, ....: ) ....: In [25]: ds.to_zarr("path/to/directory.zarr")Out[25]: <xarray.backends.zarr.ZarrStore at 0x7fe23de294c0>
(The suffix .zarr is optional–just a reminder that a zarr store livesthere.) If the directory does not exist, it will be created. If a zarrstore is already present at that path, an error will be raised, preventing itfrom being overwritten. To override this behavior and overwrite an existingstore, add mode='w' when invoking to_zarr().
DataArrays can also be saved to disk using the DataArray.to_zarr() method,and loaded from disk using the open_dataarray() function with engine=’zarr’.Similar to DataArray.to_netcdf(), DataArray.to_zarr() willconvert the DataArray to a Dataset before saving, and then convert backwhen loading, ensuring that the DataArray that is loaded is always exactlythe same as the one that was saved.
Note
xarray does not write NCZarr attributes.Therefore, NCZarr data must be opened in read-only mode.
To store variable length strings, convert them to object arrays first withdtype=object.
To read back a zarr dataset that has been created this way, we use theopen_zarr() method:
In [26]: ds_zarr = xr.open_zarr("path/to/directory.zarr")In [27]: ds_zarrOut[27]: <xarray.Dataset> Size: 248BDimensions: (x: 4, y: 5)Coordinates: * x (x) int64 32B 10 20 30 40 * y (y) datetime64[ns] 40B 2000-01-01 2000-01-02 ... 2000-01-05 z (x) <U1 16B dask.array<chunksize=(4,), meta=np.ndarray>Data variables: foo (x, y) float64 160B dask.array<chunksize=(4, 5), meta=np.ndarray>
Cloud Storage Buckets#
It is possible to read and write xarray datasets directly from / to cloudstorage buckets using zarr. This example uses the gcsfs package to providean interface to Google Cloud Storage.
General fsspec URLs, those that begin with s3:// or gcs:// for example,are parsed and the store set up for you automatically when reading.You should include any arguments to the storage backend as thekey `storage_options, part of backend_kwargs.
ds_gcs = xr.open_dataset( "gcs://<bucket-name>/path.zarr", backend_kwargs={ "storage_options": {"project": "<project-name>", "token": None} }, engine="zarr",)
This also works with open_mfdataset, allowing you to pass a list of paths ora URL to be interpreted as a glob string.
For writing, you must explicitly set up a MutableMappinginstance and pass this, as follows:
import gcsfsfs = gcsfs.GCSFileSystem(project="<project-name>", token=None)gcsmap = gcsfs.mapping.GCSMap("<bucket-name>", gcs=fs, check=True, create=False)# write to the bucketds.to_zarr(store=gcsmap)# read it backds_gcs = xr.open_zarr(gcsmap)
(or use the utility function fsspec.get_mapper()).
Zarr Compressors and Filters#
There are many different options for compression and filtering possible withzarr.
These options can be passed to the to_zarr method as variable encoding.For example:
In [28]: import zarrIn [29]: compressor = zarr.Blosc(cname="zstd", clevel=3, shuffle=2)In [30]: ds.to_zarr("foo.zarr", encoding={"foo": {"compressor": compressor}})Out[30]: <xarray.backends.zarr.ZarrStore at 0x7fe23f27bbc0>
Note
Not all native zarr compression and filtering options have been tested withxarray.
Consolidated Metadata#
Xarray needs to read all of the zarr metadata when it opens a dataset.In some storage mediums, such as with cloud object storage (e.g. Amazon S3),this can introduce significant overhead, because two separate HTTP calls to theobject store must be made for each variable in the dataset.By default Xarray uses a feature calledconsolidated metadata, storing all metadata for the entire dataset with asingle key (by default called .zmetadata). This typically drastically speedsup opening the store. (For more information on this feature, consult thezarr docs on consolidating metadata.)
By default, xarray writes consolidated metadata and attempts to read storeswith consolidated metadata, falling back to use non-consolidated metadata forreads. Because this fall-back option is so much slower, xarray issues aRuntimeWarning with guidance when reading with consolidated metadata fails:
Failed to open Zarr store with consolidated metadata, falling back to tryreading non-consolidated metadata. This is typically much slower foropening a dataset. To silence this warning, consider:
Consolidating metadata in this existing store with
zarr.consolidate_metadata().Explicitly setting
consolidated=False, to avoid trying to readconsolidate metadata.Explicitly setting
consolidated=True, to raise an error in this caseinstead of falling back to try reading non-consolidated metadata.
Modifying existing Zarr stores#
Xarray supports several ways of incrementally writing variables to a Zarrstore. These options are useful for scenarios when it is infeasible orundesirable to write your entire dataset at once.
Use
mode='a'to add or overwrite entire variables,Use
append_dimto resize and append to existing variables, andUse
regionto write to limited regions of existing arrays.
Tip
For Dataset objects containing dask arrays, asingle call to to_zarr() will write all of your data in parallel.
Warning
Alignment of coordinates is currently not checked when modifying anexisting Zarr store. It is up to the user to ensure that coordinates areconsistent.
To add or overwrite entire variables, simply call to_zarr()with mode='a' on a Dataset containing the new variables, passing in anexisting Zarr store or path to a Zarr store.
To resize and then append values along an existing dimension in a store, setappend_dim. This is a good option if data always arrives in a particularorder, e.g., for time-stepping a simulation:
In [31]: ds1 = xr.Dataset( ....: {"foo": (("x", "y", "t"), np.random.rand(4, 5, 2))}, ....: coords={ ....: "x": [10, 20, 30, 40], ....: "y": [1, 2, 3, 4, 5], ....: "t": pd.date_range("2001-01-01", periods=2), ....: }, ....: ) ....: In [32]: ds1.to_zarr("path/to/directory.zarr")Out[32]: <xarray.backends.zarr.ZarrStore at 0x7fe23f27af40>In [33]: ds2 = xr.Dataset( ....: {"foo": (("x", "y", "t"), np.random.rand(4, 5, 2))}, ....: coords={ ....: "x": [10, 20, 30, 40], ....: "y": [1, 2, 3, 4, 5], ....: "t": pd.date_range("2001-01-03", periods=2), ....: }, ....: ) ....: In [34]: ds2.to_zarr("path/to/directory.zarr", append_dim="t")Out[34]: <xarray.backends.zarr.ZarrStore at 0x7fe23f278640>
Finally, you can use region to write to limited regions of existing arraysin an existing Zarr store. This is a good option for writing data in parallelfrom independent processes.
To scale this up to writing large datasets, the first step is creating aninitial Zarr store without writing all of its array data. This can be done byfirst creating a Dataset with dummy values stored in dask,and then calling to_zarr with compute=False to write only metadata(including attrs) to Zarr:
In [35]: import dask.array# The values of this dask array are entirely irrelevant; only the dtype,# shape and chunks are usedIn [36]: dummies = dask.array.zeros(30, chunks=10)In [37]: ds = xr.Dataset({"foo": ("x", dummies)}, coords={"x": np.arange(30)})In [38]: path = "path/to/directory.zarr"# Now we write the metadata without computing any array valuesIn [39]: ds.to_zarr(path, compute=False)Out[39]: Delayed('_finalize_store-1c9cf17e-f7e4-4bf5-9ceb-f38d348f894f')
Now, a Zarr store with the correct variable shapes and attributes exists thatcan be filled out by subsequent calls to to_zarr.Setting region="auto" will open the existing store and determine thecorrect alignment of the new data with the existing coordinates, or as anexplicit mapping from dimension names to Python slice objects indicatingwhere the data should be written (in index space, not label space), e.g.,
# For convenience, we'll slice a single dataset, but in the real use-case# we would create them separately possibly even from separate processes.In [40]: ds = xr.Dataset({"foo": ("x", np.arange(30))}, coords={"x": np.arange(30)})# Any of the following region specifications are validIn [41]: ds.isel(x=slice(0, 10)).to_zarr(path, region="auto")Out[41]: <xarray.backends.zarr.ZarrStore at 0x7fe23f27b040>In [42]: ds.isel(x=slice(10, 20)).to_zarr(path, region={"x": "auto"})Out[42]: <xarray.backends.zarr.ZarrStore at 0x7fe23f27b0c0>In [43]: ds.isel(x=slice(20, 30)).to_zarr(path, region={"x": slice(20, 30)})Out[43]: <xarray.backends.zarr.ZarrStore at 0x7fe23f27b240>
Concurrent writes with region are safe as long as they modify distinctchunks in the underlying Zarr arrays (or use an appropriate lock).
As a safety check to make it harder to inadvertently override existing values,if you set region then all variables included in a Dataset must havedimensions included in region. Other variables (typically coordinates)need to be explicitly dropped and/or written in a separate calls to to_zarrwith mode='a'.
Specifying chunks in a zarr store#
Chunk sizes may be specified in one of three ways when writing to a zarr store:
Manual chunk sizing through the use of the
encodingargument in Dataset.to_zarr():Automatic chunking based on chunks in dask arrays
Default chunk behavior determined by the zarr library
The resulting chunks will be determined based on the order of the above list; daskchunks will be overridden by manually-specified chunks in the encoding argument,and the presence of either dask chunks or chunks in the encoding attribute willsupersede the default chunking heuristics in zarr.
Importantly, this logic applies to every array in the zarr store individually,including coordinate arrays. Therefore, if a dataset contains one or more daskarrays, it may still be desirable to specify a chunk size for the coordinate arrays(for example, with a chunk size of -1 to include the full coordinate).
To specify chunks manually using the encoding argument, provide a nesteddictionary with the structure {'variable_or_coord_name': {'chunks': chunks_tuple}}.
Note
The positional ordering of the chunks in the encoding argument must match thepositional ordering of the dimensions in each array. Watch out for arrays withdifferently-ordered dimensions within a single Dataset.
For example, let’s say we’re working with a dataset with dimensions('time', 'x', 'y'), a variable Tair which is chunked in x and y,and two multi-dimensional coordinates xc and yc:
In [44]: ds = xr.tutorial.open_dataset("rasm")In [45]: ds["Tair"] = ds["Tair"].chunk({"x": 100, "y": 100})In [46]: dsOut[46]: <xarray.Dataset> Size: 17MBDimensions: (time: 36, y: 205, x: 275)Coordinates: * time (time) object 288B 1980-09-16 12:00:00 ... 1983-08-17 00:00:00 xc (y, x) float64 451kB ... yc (y, x) float64 451kB ...Dimensions without coordinates: y, xData variables: Tair (time, y, x) float64 16MB dask.array<chunksize=(36, 100, 100), meta=np.ndarray>Attributes: title: /workspace/jhamman/processed/R1002RBRxaaa01a/l... institution: U.W. source: RACM R1002RBRxaaa01a output_frequency: daily output_mode: averaged convention: CF-1.4 references: Based on the initial model of Liang et al., 19... comment: Output from the Variable Infiltration Capacity... nco_openmp_thread_number: 1 NCO: netCDF Operators version 4.7.9 (Homepage = htt... history: Fri Aug 7 17:57:38 2020: ncatted -a bounds,,d...
These multi-dimensional coordinates are only two-dimensional and take up very littlespace on disk or in memory, yet when writing to disk the default zarr behavior is tosplit them into chunks:
In [47]: ds.to_zarr("path/to/directory.zarr", mode="w")Out[47]: <xarray.backends.zarr.ZarrStore at 0x7fe23f27b440>In [48]: ! ls -R path/to/directory.zarrpath/to/directory.zarr:Tair time xcycpath/to/directory.zarr/Tair:0.0.0 0.0.1 0.0.2 0.1.0 0.1.1 0.1.2 0.2.0 0.2.10.2.2path/to/directory.zarr/time:0path/to/directory.zarr/xc:0.0 1.0path/to/directory.zarr/yc:0.0 1.0
This may cause unwanted overhead on some systems, such as when reading from a cloudstorage provider. To disable this chunking, we can specify a chunk size equal to thelength of each dimension by using the shorthand chunk size -1:
In [49]: ds.to_zarr( ....: "path/to/directory.zarr", ....: encoding={"xc": {"chunks": (-1, -1)}, "yc": {"chunks": (-1, -1)}}, ....: mode="w", ....: ) ....: Out[49]: <xarray.backends.zarr.ZarrStore at 0x7fe23e443340>In [50]: ! ls -R path/to/directory.zarrpath/to/directory.zarr:Tair time xcycpath/to/directory.zarr/Tair:0.0.0 0.0.1 0.0.2 0.1.0 0.1.1 0.1.2 0.2.0 0.2.10.2.2path/to/directory.zarr/time:0path/to/directory.zarr/xc:0.0path/to/directory.zarr/yc:0.0
The number of chunks on Tair matches our dask chunks, while there is now only a singlechunk in the directory stores of each coordinate.
Iris#
The Iris tool allows easy reading of common meteorological and climate model formats(including GRIB and UK MetOffice PP files) into Cube objects which are in many ways verysimilar to DataArray objects, while enforcing a CF-compliant data model. If iris isinstalled, xarray can convert a DataArray into a Cube usingDataArray.to_iris():
In [51]: da = xr.DataArray( ....: np.random.rand(4, 5), ....: dims=["x", "y"], ....: coords=dict(x=[10, 20, 30, 40], y=pd.date_range("2000-01-01", periods=5)), ....: ) ....: In [52]: cube = da.to_iris()In [53]: cubeOut[53]: <iris 'Cube' of unknown / (unknown) (x: 4; y: 5)>
Conversely, we can create a new DataArray object from a Cube usingDataArray.from_iris():
In [54]: da_cube = xr.DataArray.from_iris(cube)In [55]: da_cubeOut[55]: <xarray.DataArray (x: 4, y: 5)> Size: 160B[20 values with dtype=float64]Coordinates: * x (x) int64 32B 10 20 30 40 * y (y) datetime64[ns] 40B 2000-01-01 2000-01-02 ... 2000-01-05
OPeNDAP#
Xarray includes support for OPeNDAP (via the netCDF4 library or Pydap), whichlets us access large datasets over HTTP.
For example, we can open a connection to GBs of weather data produced by thePRISM project, and hosted by IRI at Columbia:
In [56]: remote_data = xr.open_dataset( ....: "http://iridl.ldeo.columbia.edu/SOURCES/.OSU/.PRISM/.monthly/dods", ....: decode_times=False, ....: ) ....: In [57]: remote_dataOut[57]: <xarray.Dataset>Dimensions: (T: 1422, X: 1405, Y: 621)Coordinates: * X (X) float32 -125.0 -124.958 -124.917 -124.875 -124.833 -124.792 -124.75 ... * T (T) float32 -779.5 -778.5 -777.5 -776.5 -775.5 -774.5 -773.5 -772.5 -771.5 ... * Y (Y) float32 49.9167 49.875 49.8333 49.7917 49.75 49.7083 49.6667 49.625 ...Data variables: ppt (T, Y, X) float64 ... tdmean (T, Y, X) float64 ... tmax (T, Y, X) float64 ... tmin (T, Y, X) float64 ...Attributes: Conventions: IRIDL expires: 1375315200
Note
Like many real-world datasets, this dataset does not entirely followCF conventions. Unexpected formats will usually cause xarray’s automaticdecoding to fail. The way to work around this is to either setdecode_cf=False in open_dataset to turn off all use of CFconventions, or by only disabling the troublesome parser.In this case, we set decode_times=False because the time axis hereprovides the calendar attribute in a format that xarray does not expect(the integer 360 instead of a string like '360_day').
We can select and slice this data any number of times, and nothing is loadedover the network until we look at particular values:
In [58]: tmax = remote_data["tmax"][:500, ::3, ::3]In [59]: tmaxOut[59]: <xarray.DataArray 'tmax' (T: 500, Y: 207, X: 469)>[48541500 values with dtype=float64]Coordinates: * Y (Y) float32 49.9167 49.7917 49.6667 49.5417 49.4167 49.2917 ... * X (X) float32 -125.0 -124.875 -124.75 -124.625 -124.5 -124.375 ... * T (T) float32 -779.5 -778.5 -777.5 -776.5 -775.5 -774.5 -773.5 ...Attributes: pointwidth: 120 standard_name: air_temperature units: Celsius_scale expires: 1443657600# the data is downloaded automatically when we make the plotIn [60]: tmax[0].plot()
Some servers require authentication before we can access the data. For thispurpose we can explicitly create a backends.PydapDataStoreand pass in a Requests session object. For example forHTTP Basic authentication:
import xarray as xrimport requestssession = requests.Session()session.auth = ('username', 'password')store = xr.backends.PydapDataStore.open('http://example.com/data', session=session)ds = xr.open_dataset(store)
Pydap’s cas module has functions that generate custom sessions forservers that use CAS single sign-on. For example, to connect to serversthat require NASA’s URS authentication:
import xarray as xrfrom pydata.cas.urs import setup_sessionds_url = 'https://gpm1.gesdisc.eosdis.nasa.gov/opendap/hyrax/example.nc'session = setup_session('username', 'password', check_url=ds_url)store = xr.backends.PydapDataStore.open(ds_url, session=session)ds = xr.open_dataset(store)
Pickle#
The simplest way to serialize an xarray object is to use Python’s built-in picklemodule:
In [61]: import pickle# use the highest protocol (-1) because it is way faster than the default# text based pickle formatIn [62]: pkl = pickle.dumps(ds, protocol=-1)In [63]: pickle.loads(pkl)Out[63]: <xarray.Dataset> Size: 17MBDimensions: (time: 36, y: 205, x: 275)Coordinates: * time (time) object 288B 1980-09-16 12:00:00 ... 1983-08-17 00:00:00 xc (y, x) float64 451kB ... yc (y, x) float64 451kB ...Dimensions without coordinates: y, xData variables: Tair (time, y, x) float64 16MB dask.array<chunksize=(36, 100, 100), meta=np.ndarray>Attributes: title: /workspace/jhamman/processed/R1002RBRxaaa01a/l... institution: U.W. source: RACM R1002RBRxaaa01a output_frequency: daily output_mode: averaged convention: CF-1.4 references: Based on the initial model of Liang et al., 19... comment: Output from the Variable Infiltration Capacity... nco_openmp_thread_number: 1 NCO: netCDF Operators version 4.7.9 (Homepage = htt... history: Fri Aug 7 17:57:38 2020: ncatted -a bounds,,d...
Pickling is important because it doesn’t require any external librariesand lets you use xarray objects with Python modules likemultiprocessing or Dask. However, pickling isnot recommended for long-term storage.
Restoring a pickle requires that the internal structure of the types for thepickled data remain unchanged. Because the internal design of xarray is stillbeing refined, we make no guarantees (at this point) that objects pickled withthis version of xarray will work in future versions.
Note
When pickling an object opened from a NetCDF file, the pickle file willcontain a reference to the file on disk. If you want to store the actualarray values, load it into memory first with Dataset.load()or Dataset.compute().
Dictionary#
We can convert a Dataset (or a DataArray) to a dict usingDataset.to_dict():
In [64]: ds = xr.Dataset({"foo": ("x", np.arange(30))})In [65]: dsOut[65]: <xarray.Dataset> Size: 240BDimensions: (x: 30)Dimensions without coordinates: xData variables: foo (x) int64 240B 0 1 2 3 4 5 6 7 8 9 ... 21 22 23 24 25 26 27 28 29In [66]: d = ds.to_dict()In [67]: dOut[67]: {'coords': {}, 'attrs': {}, 'dims': {'x': 30}, 'data_vars': {'foo': {'dims': ('x',), 'attrs': {}, 'data': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]}}}
We can create a new xarray object from a dict usingDataset.from_dict():
In [68]: ds_dict = xr.Dataset.from_dict(d)In [69]: ds_dictOut[69]: <xarray.Dataset> Size: 240BDimensions: (x: 30)Dimensions without coordinates: xData variables: foo (x) int64 240B 0 1 2 3 4 5 6 7 8 9 ... 21 22 23 24 25 26 27 28 29
Dictionary support allows for flexible use of xarray objects. It doesn’trequire external libraries and dicts can easily be pickled, or converted tojson, or geojson. All the values are converted to lists, so dicts mightbe quite large.
To export just the dataset schema without the data itself, use thedata=False option:
In [70]: ds.to_dict(data=False)Out[70]: {'coords': {}, 'attrs': {}, 'dims': {'x': 30}, 'data_vars': {'foo': {'dims': ('x',), 'attrs': {}, 'dtype': 'int64', 'shape': (30,)}}}
This can be useful for generating indices of dataset contents to expose tosearch indices or other automated data discovery tools.
Rasterio#
GDAL readable raster data using rasterio such as GeoTIFFs can be opened using the rioxarray extension.rioxarray can also handle geospatial related tasks such as re-projecting and clipping.
In [71]: import rioxarrayIn [72]: rds = rioxarray.open_rasterio("RGB.byte.tif")In [73]: rdsOut[73]: <xarray.DataArray (band: 3, y: 718, x: 791)>[1703814 values with dtype=uint8]Coordinates: * band (band) int64 1 2 3 * y (y) float64 2.827e+06 2.826e+06 ... 2.612e+06 2.612e+06 * x (x) float64 1.021e+05 1.024e+05 ... 3.389e+05 3.392e+05 spatial_ref int64 0Attributes: STATISTICS_MAXIMUM: 255 STATISTICS_MEAN: 29.947726688477 STATISTICS_MINIMUM: 0 STATISTICS_STDDEV: 52.340921626611 transform: (300.0379266750948, 0.0, 101985.0, 0.0, -300.0417827... _FillValue: 0.0 scale_factor: 1.0 add_offset: 0.0 grid_mapping: spatial_refIn [74]: rds.rio.crsOut[74]: CRS.from_epsg(32618)In [75]: rds4326 = rds.rio.reproject("epsg:4326")In [76]: rds4326.rio.crsOut[76]: CRS.from_epsg(4326)In [77]: rds4326.rio.to_raster("RGB.byte.4326.tif")
GRIB format via cfgrib#
Xarray supports reading GRIB files via ECMWF cfgrib python driver,if it is installed. To open a GRIB file supply engine='cfgrib'to open_dataset() after installing cfgrib:
In [78]: ds_grib = xr.open_dataset("example.grib", engine="cfgrib")
We recommend installing cfgrib via conda:
conda install -c conda-forge cfgrib
CSV and other formats supported by pandas#
For more options (tabular formats and CSV files in particular), considerexporting your objects to pandas and using its broad range of IO tools.For CSV files, one might also consider xarray_extras.
Third party libraries#
More formats are supported by extension libraries:
xarray-mongodb: Store xarray objects on MongoDB
