Reviving this exploration with xarray : 2023.3.0 and I'm no longer seeing the traceback linked above, so I think this can be closed

Original traceback from xr.align(copy, newds), which now succeeds: ``pytbFile ~/mambaforge/envs/xarray-release/lib/python3.10/site-packages/xarray/core/indexes.py:490, in PandasIndex.equals(self, other) 488 if not isinstance(other, PandasIndex): 489 return False --> 490 return self.index.equals(other.index) and self.dim == other.dim

AttributeError: 'PandasIndex' object has no attribute 'index' ```

@leomiquelutti I wasn't able to reproduce this (ran your colab notebook and opened the netcdf in qgis 3.26). It's possible for versions to change between your colab executions so I recommend using rioxarray.show_versions() Since you are using rioxarray to write the CRS, I'd also recommend opening this issue or discussion over in https://github.com/corteva/rioxarray

This would be a pretty small change and only applies for loading data into numpy arrays, for example current repr for a variable followed by modified for the example dataset above (which already happens for large arrays):

Seeing a few values at the edges can be nice, so this makes me realize how data summaries in the metadata (Zarr or STAC) is great for large datasets on cloud storage.

One of the main motivations behind the the rioxarray extension is GDAL compatibility. It looks like @snowman2 and @TomAugspurger have discussed saving many geotiffs loaded into xarray as GDAL-compatible Zarr for example https://github.com/corteva/rioxarray/issues/433#issuecomment-967685356.

While it seems that the ultimate solution is agreeing on a format standard, here is another small example using the rioxarray extension where format conversion doesn't currently work as you might expect:

```python

https://github.com/pydata/xarray-data

ds = xr.open_dataset('xarray-data/air_temperature.nc', engine='rasterio')

TooManyDimensions: Only 2D and 3D data arrays supported.

ds.rio.to_raster('test.zarr', driver='ZARR')

Does not error, but output not equivalent to gdal_translate -of ZARR xarray-data/air_temperature.nc gdal_air_temp.zarr

for example, gdalinfo xarray-tutorial-airtemp.zarr gives

Warning 1: Too many samples along the > 2D dimensions of /air.

ds.to_zarr('xarray-tutorial-airtemp.zarr')

```

Well, I changed that URL back in https://github.com/pydata/xarray/pull/3162 . But maybe it's best to just remove that test now that VRT and other GDAL/rasterio functionality is integrated with rioxarray?

Unfortunately the AWS "landsat-pds" (public data set) was deprecated and not supported any longer, so maybe those images are finally being deleted https://lists.osgeo.org/pipermail/landsat-pds/2021-July/000185.html

Thanks @keewis I should have been more clear about the environment. I was recently going over the tutorial with someone and started with:

1. conda create -n xarray-tutorial xarray running ds = xr.tutorial.load_dataset("air_temperature") --> ImportError: using the tutorial data requires pooch
2. we install pooch and then hit: ValueError: cannot guess the engine, try passing one explicitly
3. after consulting the docstring we then try ds = xr.tutorial.load_dataset("air_temperature", engine="netcdf4") and hit ValueError: unrecognized engine netcdf4 must be one of: ['store']
4. being familiar with xarray we then installed netcdf4 into our environment and all is well.

I do think these error messages are not obvious to fix for new xarray users trying out the tutorial (especially # 3 above)

we could definitely improve the error message, though. Something like "unknown engine {engine}, please choose one of the installed engines: {engines}", maybe?

Yes. Perhaps with a link to docs with a list of engines? for the tutorial case specifically could also update the ImportError message to read ImportError: please install 'pooch' and 'netcdf4' to use xarray tutorial data?

Just wanted to rekindle discussion here and ping @dcherian and @benbovy , the current workaround for pandas DatetimeIndex with timezone info (dtype='datetime64[ns, EST]') is to drop the timezone piece or use to_index() and operate in pandas, then reassign the time coordinate: See https://github.com/pydata/xarray/issues/1036 and https://github.com/pydata/xarray/issues/3163.

If I'm following https://github.com/pydata/xarray/blob/master/design_notes/flexible_indexes_notes.md this is another potential example of improved user-friendliness where we could have timezone-aware indexes and therefore call pandas methods like pandas.core.indexes.datetimes.DatetimeIndex.tz_convert() directly as a DataArray method?

This would definitely be great for remote sensing data that is usually stored with UTC timestamps, but often analysis requires converting to local time.

Thanks @gabriel-abrahao good catch! the docs are definitely outdated here and need to be changed. This has been brought up before here https://github.com/pydata/xarray/issues/3185#issuecomment-574215734

You can either do what you suggested or simply transform = Affine(*da.transform) https://github.com/pydata/xarray/blob/12b4480ff2bde696142ca850275cdcc85ca0fbc9/xarray/tests/test_backends.py#L4468

A PR to update the documentation would be welcome. It'd be just a one-line fix here: https://github.com/pydata/xarray/blob/d2582c2f8811a3bd527d47c945b1cccd4983a1d3/xarray/backends/rasterio_.py#L179

cc @snowman2 the docstring for open_rasterio has the same issue currently in rioxarray https://github.com/corteva/rioxarray/blob/63ce87f92835897d46aaa270f5ec07aca4e72cda/rioxarray/_io.py#L653

@Tinkaa I looked into this a bit more and I suspect the way you are installing packages (conda, pip?) is important to bring in compatible libnetcdf behind the scenes (You're xr.show_versions() shows new libnetcdf: 4.7.4 but older netCDF4: 1.5.0.1 making me think you installed and older version netcdf into an existing environment without any dependency resolution).

This works for me on mac0S: conda create -c conda-forge -n opendap_working python xarray netcdf4==1.5.1 matplotlib Note that I found the explicit pin important. Perhaps because it brings in libnetcdf=4.6.2 behind the scenes? Here is my full list of working versions:

``` INSTALLED VERSIONS

commit: None python: 3.7.10 | packaged by conda-forge | (default, Feb 19 2021, 15:59:12) [Clang 11.0.1 ] python-bits: 64 OS: Darwin OS-release: 20.3.0 machine: x86_64 processor: i386 byteorder: little LC_ALL: None LANG: en_US.UTF-8 LOCALE: en_US.UTF-8 libhdf5: 1.10.4 libnetcdf: 4.6.2

xarray: 0.16.2 pandas: 1.2.2 numpy: 1.20.1 scipy: None netCDF4: 1.5.0.1 pydap: None h5netcdf: None h5py: None Nio: None zarr: None cftime: 1.4.1 nc_time_axis: None PseudoNetCDF: None rasterio: None cfgrib: None iris: None bottleneck: None dask: None distributed: None matplotlib: 3.3.4 cartopy: None seaborn: None numbagg: None pint: None setuptools: 49.6.0.post20210108 pip: 21.0.1 conda: None pytest: None IPython: None sphinx: None ```

Hopefully some folks that know more about opendap could chime in (@dopplershift, @rabernat ?) to help get to the root of this issue. For what it's worth, I noticed this additional log output that wasn't surfaced in the jupyter notebook in addition to the original traceback i posted: oc_open: server error retrieving url: code=3 message="The identifier `tmax.tmax%5b0%5d%5b0:3:620%5d%5b0:3:1404%5d' is not in the dataset."oc_open: server error retrieving url: code=3 message="The identifier `tmax.tmax%5b0%5d%5b0:3:620%5d%5b0:3:1404%5d' is not in the dataset."oc_open: server error retrieving url: code=3 message="The identifier `tmax.tmax%5b0%5d%5b0:3:620%5d%5b0:3:1404%5d' is not in the dataset."oc_open: server error retrieving url: code=3 message="The identifier `tmax.tmax%5b0%5d%5b0:3:620%5d%5b0:3:1404%5d' is not in the dataset." oc_open: server error retrieving url: code=3 message="The identifier `tmax.tmax%5b0%5d%5b0:3:620%5d%5b0:3:1404%5d' is not in the dataset."^CTraceback (most recent call last): File "/Users/scott/miniconda3/envs/opendap_test/lib/python3.9/site-packages/xarray/backends/common.py", line 52, in robust_getitem return array[key] File "src/netCDF4/_netCDF4.pyx", line 4420, in netCDF4._netCDF4.Variable.__getitem__ File "src/netCDF4/_netCDF4.pyx", line 5363, in netCDF4._netCDF4.Variable._get File "src/netCDF4/_netCDF4.pyx", line 1950, in netCDF4._netCDF4._ensure_nc_success RuntimeError: NetCDF: file not found

thanks for the ping @dcherian, i really like the idea! One other thing that often gets neglected in test suites is operating on remote data. I understand the need to avoid long-running tests and tests prone to network failures for PRs, but running these sorts of examples as a cron job could be very helpful for benchmarking and detecting issues.

In intake-xarray we recently added tests against a local HTTP server and "S3" server: https://github.com/intake/intake-xarray/blob/master/intake_xarray/tests/test_remote.py

Also added several simple tests requiring a network connection to public data (no auth required) that we run locally but not in CI currently: https://github.com/intake/intake-xarray/blob/master/intake_xarray/tests/test_network.py

@dcherian I was just revisiting a use-case for this and realized xr.open_rasterio() does not take **kwargs. http://xarray.pydata.org/en/stable/generated/xarray.open_rasterio.html#xarray.open_rasterio

I think this would be easy to implement because a while back rasterio implemented a keyword argument so that you can do rasterio.open(path, overview_level=3) https://github.com/mapbox/rasterio/issues/1504

So is this just a matter of accepting kwargs in xr.open_rasterio() and passing them through? https://github.com/pydata/xarray/blob/4aa7622b6ff16647df64fe69f39438b7cbe9576c/xarray/backends/rasterio_.py#L241

As seems to be done here for example: https://github.com/pydata/xarray/blob/cc271e61077c543e0f3b1a06ad5e905ea2c91617/xarray/backends/h5netcdf_.py#L164

Not sure, but I think the h5netcdf engine is the only one that allows for file-like objects (so anything going through fsspec)

Awesome @dcherian ! I like the approach of defining a binder environment then pulling the examples directory into the binder session with nbgitpuller. There are options for where to store the binder environment config (1) master branch .binder, 2) a new branch called 'binder', or 3) even a separate repo... maybe pydata/xarray-binder-env - see https://github.com/scottyhq/repo2docker-githubci ?

ok. as noted https://github.com/mapbox/rasterio/pull/1709 the aws_unsigned=True should bypass importing boto3 in upcoming rasterio 1.0.25. But I've gone ahead and added boto3 to ci requirements.yml files in case people want to add other tests reading from s3, gcs, etc. in the future

I'm in no rush to merge this (just wanted to help get rid of the red on Azure :). I think it's important to have network tests since it's increasingly common to be streaming data from s3, gcs, etc rather than reading locally. But I'm not sure of best practices for setting such tests up, so any suggestions for modifications are welcome.

but seems like the tests are now run with each build on Azure? https://dev.azure.com/xarray/xarray/_build/results?buildId=389&view=ms.vss-test-web.build-test-results-tab

the new test currently still fails - I think due to https://github.com/mapbox/rasterio/pull/1709

cc @shoyer who created the issue linked in the first post.

sounds good. i'll try to get to it this week.

Strange, the same test works on my laptop current packages. Some ideas below:

• in the rasterio.Envenvironment we can add CPL_CURL_VERBOSE=True to get more information about the error, but would need to capture output that typically goes to the terminal (I've done this before with https://github.com/minrk/wurlitzer).

• I'm inclined to remove CPL_VSIL_CURL_USE_HEAD=False which is not documented anywhere obvious (https://lists.osgeo.org/pipermail/gdal-dev/2014-August/039924.html) and could be behind the error message. All rasterio/gdal environment config options are detailed here https://trac.osgeo.org/gdal/wiki/ConfigOptions

• we might also consider a different URL. maybe add a tif file to https://github.com/pydata/xarray-data? or use one of same URLs that the rasterio library is using for network tif tests? https://github.com/mapbox/rasterio/blob/c2a495d536db6328e018cb65ce7cc0b8c559a937/tests/test_env.py#L34-L36

@jmichel-otb completely agree on your comments and think that the remote sensing community gains a lot by making these packages work well together (gdal, rasterio, xarray, dask, etc). As more people start using the rasterio backend with various use cases hopefully there will be more contributions such as this. Thanks!

thanks for the input @shoyer, I attempted to tidy up a bit and in the process re-ordered some things such as adding an 'engine' check at the top of open_dataset(). backend tests are passing locally on my machine. hopefully i didn't add too much here or overstep!

@shoyer , it would be great to have your feedback on these recent changes now that h5netcdf 0.7 is out. There's a bit more logic required in api.py now that scipy isn't the only backend that is able to read file-like objects (and people may not specify engine= when opening datasets)

test_backends.py passes locally for me except for TestValidateAttrs.test_validating_attrs... not sure why.

Also, per your comment here: https://github.com/shoyer/h5netcdf/pull/51#issuecomment-467591446, I think it would be great to get a few small netcdf4/hdf test files in https://github.com/pydata/xarray-data.

Just noting here that I've gotten this to work reading a netcd4/fhdf5 file via gcsfs, but not for the same file accessed via s3fs: https://github.com/dask/s3fs/issues/168

I think the minimum rasterio version should be increased to 1.0. For whatever reason the conda defaults channel hasn't been updated since 0.36 (Jun 14, 2016!). There are many important changes in 1.0 and beyond, and those releases are available via both pip and conda-forge.

thanks for the feedback @fmaussion, I think I've addressed your suggestions, let me know if anything else needs adjusting

Following up on https://github.com/dask/dask/issues/3255 @mrocklin, @shoyer, @jhamman

As @darothen mentioned, first thing is to check that the geotiffs themselves are tiled (otherwise I'm guessing that open_rasterio() will open the entire thing. You can do this with:

python import rasterio with rasterio.open('image_001.tif') as src: print(src.profile)

Here is the mentioned example notebook which works for tiled geotiffs stored on google cloud: https://github.com/scottyhq/pangeo-example-notebooks/tree/binderfy

You can use the 'launch binder' button to run it with a pangeo dask-kubernetes cluster, or just read through the landsat8-cog-ndvi.ipynb notebook.