issue_comments

the join allow some float imprecision (similar to method=nearest), which would conveniently allow cases like this to work

I like that.

@shoyer, I'd volunteer for a PR, should you agree extending Dataset.to_zarr in a backward compatible way:

python def to_zarr(self, ..., encode_cf: bool = True): ...

After debugging we found that zarr.core.Array._encode_chunk() does not encode chunks, if both compressor and filters are missing,

However I could not reproduce our problem with Zarr open/save alone. It seems to occur only when using xarrays.open_zarr() and xr.Dataset.to_zarr(). Therefore I seems to be an xarray issue rather than a Zarr one.

I'm still suffering from IndexError: pop from an empty deque. Can somebody tell me which s3fs version to use after fix by @martindurant?

Here are my relevant packages:

# Name                    Version                   Build  Channel
aiobotocore               0.10.3                     py_0    conda-forge
aiohttp                   3.7.3            py39hb82d6ee_0    conda-forge
botocore                  1.12.91                    py_0    conda-forge
dask                      2.30.0                     py_0    conda-forge
dask-core                 2.30.0                     py_0    conda-forge
distributed               2.30.1           py39hcbf5309_0    conda-forge
fsspec                    0.8.4                      py_0    conda-forge
python                    3.9.0           h7840368_5_cpython    conda-forge
s3fs                      0.5.1                      py_0    conda-forge
xarray                    0.16.2             pyhd8ed1ab_0    conda-forge
zarr                      2.5.0                      py_0    conda-forge

Thanks in advance!

@hans-permana your example shows a different issue: indexing with a date string yields a time dimension of length 1, rather than squeezing it out

Nope, look at the screenshot again, the dimension is zero. The very similar issue (if not same) remains and should be considered a bug:

If I now use sel() with a date string without time component, I get a 3D array with zero time dimension:

However, if I use sel() with a date string with time component, I get the expected 2D array:

EDIT

It seems that if I create the cube dataset from above with a time coordinate variable whose values don't have a time component (e.g. 2018-06-26 00:00:00.000000), then both sel(time='2018-06-26') and sel(time='2018-06-26 10:23:05') work as expected and only yield 2D results.

EDIT 2

Root cause may be related to Pandas indexing using strings that encode different accuracy / resolution: http://pandas-docs.github.io/pandas-docs-travis/user_guide/timeseries.html#slice-vs-exact-match. Very contra-intuitive.

Output of xr.show_versions()

@shoyer thanks, time of the call dropped down to a second!

Yes, so this is a duplicate of https://github.com/pydata/xarray/issues/1444, sorry!

When can we expect 0.9.7 with the fix?

@shoyer

We currently decode anything with a _FillValue attribute to float, ...

I believe this fact is surprising for any user of integer/index/enum/classification datasets. Since its justification seems to be an implementation detail which comes at the cost of increased memory and CPU consumption I suggest documenting it in open_dataset() and decode_cf() functions.

Here is how we overcome this issue by deleting the _FillValue attribute of integer variables if their scale_factor and add_offset attributes are not provided:

ds = xr.open_dataset(path, decode_cf=False)
old_fill_values = unset_fill_value_for_int_vars(ds)
ds = xr.decode_cf(ds)
reset_fill_value_for_int_vars(ds, old_fill_values)

where old_fill_values is a mapping of variable names to fill values.

@jhamman _NoFill is about optimizing writes, see nc_set_fill

I see, that is what is done in mask_and_scale(). Why can't xarray used masked arrays, that would retain the original dtype? (Dask, I guess?) Expanding integers to 8 byte floats not only cost memory but also CPU, including an inaccurate in-memory integer representation.

I guess, the poblem is caused in xarray/conventions.py.

Note, when debugging into it, fill_value == nd.array([0], dtype == np.int8) and fill_value.dtype.kind='i' and the latter kind is not dealt with. Therefore int8 is turned into float64.

Here you are

$ ncdump -h -s
netcdf ESACCI-LC-L4-LCCS-Map-300m-P5Y-2005-v1.6.1 {
dimensions:
        lat = 64800 ;
        lon = 129600 ;
variables:
        byte lccs_class(lat, lon) ;
                lccs_class:long_name = "Land cover class defined in LCCS" ;
                lccs_class:standard_name = "land_cover_lccs" ;
                lccs_class:flag_values = 0b, 10b, 11b, 12b, 20b, 30b, 40b, 50b, 60b, 61b, 62b, 70b, 71b, 72b, 80b, 81b, 82b, 90b, 100b, 110b, 120b, 121b, 122b, -126b, -116b, -106b, -104b, -103b, -96b, -86b, -76b, -66b, -56b, -55b, -54b, -46b,
-36b ;
                lccs_class:flag_meanings = "no_data cropland_rainfed cropland_rainfed_herbaceous_cover cropland_rainfed_tree_or_shrub_cover cropland_irrigated mosaic_cropland mosaic_natural_vegetation tree_broadleaved_evergreen_closed_to_open
tree_broadleaved_deciduous_closed_to_open tree_broadleaved_deciduous_closed tree_broadleaved_deciduous_open tree_needleleaved_evergreen_closed_to_open tree_needleleaved_evergreen_closed tree_needleleaved_evergreen_open tree_needleleaved_decidu
ous_closed_to_open tree_needleleaved_deciduous_closed tree_needleleaved_deciduous_open tree_mixed mosaic_tree_and_shrub mosaic_herbaceous shrubland shrubland_evergreen shrubland_deciduous grassland lichens_and_mosses sparse_vegetation sparse_s
hrub sparse_herbaceous tree_cover_flooded_fresh_or_brakish_water tree_cover_flooded_saline_water shrub_or_herbaceous_cover_flooded urban bare_areas bare_areas_consolidated bare_areas_unconsolidated water snow_and_ice" ;
                lccs_class:valid_min = 1 ;
                lccs_class:valid_max = 220 ;
                lccs_class:_Unsigned = "true" ;
                lccs_class:_FillValue = 0b ;
                lccs_class:ancillary_variables = "processed_flag current_pixel_state observation_count algorithmic_confidence_level" ;
                lccs_class:_Storage = "chunked" ;
                lccs_class:_ChunkSizes = 2048, 2048 ;
                lccs_class:_DeflateLevel = 6 ;
                lccs_class:_NoFill = "true" ;
        byte processed_flag(lat, lon) ;
                processed_flag:standard_name = "land_cover_lccs status_flag" ;
                processed_flag:flag_values = 0b, 1b ;
                processed_flag:flag_meanings = "not_processed processed" ;
                processed_flag:valid_min = 0 ;
                processed_flag:valid_max = 1 ;
                processed_flag:_FillValue = -1b ;
                processed_flag:long_name = "LC map processed area flag" ;
                processed_flag:_Storage = "chunked" ;
                processed_flag:_ChunkSizes = 2048, 2048 ;
                processed_flag:_DeflateLevel = 6 ;
                processed_flag:_NoFill = "true" ;
        byte current_pixel_state(lat, lon) ;
                current_pixel_state:standard_name = "land_cover_lccs status_flag" ;
                current_pixel_state:flag_values = 0b, 1b, 2b, 3b, 4b, 5b ;
                current_pixel_state:flag_meanings = "invalid clear_land clear_water clear_snow_ice cloud cloud_shadow" ;
                current_pixel_state:valid_min = 0 ;
                current_pixel_state:valid_max = 5 ;
                current_pixel_state:_FillValue = -1b ;
                current_pixel_state:long_name = "LC pixel type mask" ;
                current_pixel_state:_Storage = "chunked" ;
                current_pixel_state:_ChunkSizes = 2048, 2048 ;
                current_pixel_state:_DeflateLevel = 6 ;
                current_pixel_state:_NoFill = "true" ;
        short observation_count(lat, lon) ;
                observation_count:standard_name = "land_cover_lccs number_of_observations" ;
                observation_count:valid_min = 0 ;
                observation_count:valid_max = 32767 ;
                observation_count:_FillValue = -1s ;
                observation_count:long_name = "number of valid observations" ;
                observation_count:_Storage = "chunked" ;
                observation_count:_ChunkSizes = 2048, 2048 ;
                observation_count:_DeflateLevel = 6 ;
                observation_count:_Endianness = "little" ;
                observation_count:_NoFill = "true" ;
        byte algorithmic_confidence_level(lat, lon) ;
                algorithmic_confidence_level:standard_name = "land_cover_lccs algorithmic_confidence" ;
                algorithmic_confidence_level:valid_min = 0 ;
                algorithmic_confidence_level:valid_max = 100 ;
                algorithmic_confidence_level:scale_factor = 0.01f ;
                algorithmic_confidence_level:_FillValue = -1b ;
                algorithmic_confidence_level:long_name = "LC map confidence level based on algorithm performance" ;
                algorithmic_confidence_level:_Storage = "chunked" ;
                algorithmic_confidence_level:_ChunkSizes = 2048, 2048 ;
                algorithmic_confidence_level:_DeflateLevel = 6 ;
                algorithmic_confidence_level:_NoFill = "true" ;
        float lat(lat) ;
                lat:long_name = "latitude" ;
                lat:standard_name = "latitude" ;
                lat:valid_min = -89.9986f ;
                lat:valid_max = 89.99861f ;
                lat:units = "degrees_north" ;
                lat:_Storage = "chunked" ;
                lat:_ChunkSizes = 64800 ;
                lat:_DeflateLevel = 6 ;
                lat:_Endianness = "little" ;
                lat:_NoFill = "true" ;
        float lon(lon) ;
                lon:long_name = "longitude" ;
                lon:standard_name = "longitude" ;
                lon:valid_min = -179.9986f ;
                lon:valid_max = 179.9986f ;
                lon:units = "degrees_east" ;
                lon:_Storage = "chunked" ;
                lon:_ChunkSizes = 129600 ;
                lon:_DeflateLevel = 6 ;
                lon:_Endianness = "little" ;
                lon:_NoFill = "true" ;
        int crs ;
                crs:i2m = "0.002777777701187,0.0,0.0,-0.002777777701187,-180.00000033927267,90.0" ;
                crs:wkt = "GEOGCS[\"WGS 84\", \r\n  DATUM[\"World Geodetic System 1984\", \r\n    SPHEROID[\"WGS 84\", 6378137.0, 298.257223563, AUTHORITY[\"EPSG\",\"7030\"]], \r\n    AUTHORITY[\"EPSG\",\"6326\"]], \r\n  PRIMEM[\"Greenwich\",
0.0, AUTHORITY[\"EPSG\",\"8901\"]], \r\n  UNIT[\"degree\", 0.017453292519943295], \r\n  AXIS[\"Geodetic longitude\", EAST], \r\n  AXIS[\"Geodetic latitude\", NORTH], \r\n  AUTHORITY[\"EPSG\",\"4326\"]]" ;
                crs:_Endianness = "little" ;
                crs:_NoFill = "true" ;

// global attributes:
                :title = "ESA CCI Land Cover Map" ;
                :summary = "This dataset contains the global ESA CCI land cover classification map derived from satellite data of one epoch." ;
                :type = "ESACCI-LC-L4-LCCS-Map-300m-P5Y" ;
                :id = "ESACCI-LC-L4-LCCS-Map-300m-P5Y-2005-v1.6.1" ;
                :project = "Climate Change Initiative - European Space Agency" ;
                :references = "http://www.esa-landcover-cci.org/" ;
                :institution = "Universite catholique de Louvain" ;
                :contact = "landcover-cci@uclouvain.be" ;
                :comment = "" ;
                :Conventions = "CF-1.6" ;
                :standard_name_vocabulary = "NetCDF Climate and Forecast (CF) Standard Names version 21" ;
                :keywords = "land cover classification,satellite,observation" ;
                :keywords_vocabulary = "NASA Global Change Master Directory (GCMD) Science Keywords" ;
                :license = "ESA CCI Data Policy: free and open access" ;
                :naming_authority = "org.esa-cci" ;
                :cdm_data_type = "grid" ;
                :TileSize = "2048:2048" ;
                :tracking_id = "00f7e0ee-3b0e-4ea3-9b9f-186e02fb4439" ;
                :product_version = "1.6.1" ;
                :date_created = "20151217T094622Z" ;
                :creator_name = "University catholique de Louvain" ;
                :creator_url = "http://www.uclouvain.be/" ;
                :creator_email = "landcover-cci@uclouvain.be" ;
                :source = "MERIS FR L1B version 5.05, MERIS RR L1B version 8.0, SPOT VGT P" ;
                :history = "amorgos-4,0, lc-sdr-1.0, lc-sr-1.0, lc-classification-1.0,lc-user-tools-3.10" ;
                :time_coverage_start = "20030101" ;
                :time_coverage_end = "20071231" ;
                :time_coverage_duration = "P5Y" ;
                :time_coverage_resolution = "P5Y" ;
                :geospatial_lat_min = "-89.99999" ;
                :geospatial_lat_max = "90.0" ;
                :geospatial_lon_min = "-180.0" ;
                :geospatial_lon_max = "179.99998" ;
                :spatial_resolution = "300m" ;
                :geospatial_lat_units = "degrees_north" ;
                :geospatial_lat_resolution = "0.002778" ;
                :geospatial_lon_units = "degrees_east" ;
                :geospatial_lon_resolution = "0.002778" ;
                :_SuperblockVersion = 2 ;
                :_IsNetcdf4 = 1 ;
                :_Format = "netCDF-4" ;
}

@PeterDSteinberg please have a look at module gridtools.resampling of repo https://github.com/CAB-LAB/gridtools. There are various up- and downsampling methods, which can deal with NaNs, and which are fast as C thanks to JIT through Numba. We use this package successfully in two projects.

I'd vote for having two functions but still have an option in xarray.concat that allows for stacking of variables by adding the new concat dimension. Same option should be available for xarray.open_mfdataset.

like we could use the bounds attribute in this dataset

Yes. And use of the bounds attribute is also CF-compliant.

Now sorry for the delay on my side - just returned from Holidays.

Here is the concrete example: https://www.dropbox.com/sh/1a30p6aya96nftl/AAD6E4aCRkC2PLafZDboFszJa?dl=0 (The *.nc files contain time series of images of analysed sea surface temperatures and are generated by the ESA SST CCI (Climate Change Initiative) project.)

If I open these using open_mfdataset() then lat_bnds, lon_bnds, time_bnds have an extra time dimension, which of course doesn't make sense.

Ideally, lat_bnds, lon_bnds, time_bnds should be correctly recognized as kind of coordinates, as you say.

Fantastic, thanks!

With h5py the data is read correctly too.

Ok, I'll submit a netCDF issue then.

Just found that the valid_min/valid_max attributes are not directly part of the CF convention but the NUG convention and applying to byte only according to CF Section 2.2 Data Types

All integer types are treated by the netCDF interface as signed. It is possible to treat the byte type as unsigned by using the NUG convention of indicating the unsigned range using the valid_min, valid_max, or valid_range attributes.

As for for #821, Panoply shows the correct values for the same file:

After some testing I found the problem. The single value of the time coordinate is wrong in the files. So it is a file content problem not a problem in the software. Therefore I'll close this issue.

However, Panoply displays the time information correctly and I found out why: Panoply correctly interprets the time_bnds variable to which the time coordinate variable points to via its attribute bounds. Since this is conforming to CF, I wonder whether xarray should support this bounds convention. From the CF conventions 1.6:

It is often the case that data values are not representative of single points in time and/or space, but rather of intervals or multidimensional cells. This convention defines a bounds attribute to specify the extent of intervals or cells. When data that is representative of cells can be described by simple statistical methods, those methods can be indicated using the cell_methods attribute. An important application of this attribute is to describe climatological and diurnal statistics.

Details are in Section 7.1 Cell Boundaries

Thanks, I'll give it a try.

@jhamman: any progress on this? Our team would be happy to contribute as we have similar requirements in our project.

Thanks for the prompt reply!

Once we have decided to use xarray for our project(s) and once we familiarized with its internals, we'll be happy to contribute and support you! Currently we all feel a bit dizzy about the many options we have and how to decide which way to go: Create our own library using xarray or build on UK MetOffice's Iris, Apache OCW, or Max-Planck-Institute's CDO, etc.