Repository: corteva/rioxarray

Issue:
`open_rasterio` support subdatasets
It would be a nice feature to be able to load in subdatasets from rasterio.

Discussion:

Patch:
diff --git a/rioxarray/_io.py b/rioxarray/_io.py
index 55f32887..8974ca49 100644
--- a/rioxarray/_io.py
+++ b/rioxarray/_io.py
@@ -12,7 +12,7 @@ from collections import OrderedDict
 from distutils.version import LooseVersion
 
 import numpy as np
-from xarray import DataArray
+from xarray import DataArray, Dataset
 from xarray.core import indexing
 from xarray.core.utils import is_scalar
 from xarray.backends.common import BackendArray
@@ -48,7 +48,7 @@ class RasterioArrayWrapper(BackendArray):
         if not np.all(np.asarray(dtypes) == dtypes[0]):
             raise ValueError("All bands should have the same dtype")
 
-        self._dtype = np.float64 if self.masked else np.dtype(dtypes[0])
+        self._dtype = np.dtype("float64") if self.masked else np.dtype(dtypes[0])
 
     @property
     def dtype(self):
@@ -174,6 +174,27 @@ def _parse_envi(meta):
     return parsed_meta
 
 
+def _parse_tags(tags):
+    def parsevec(s):
+        return np.fromstring(s.strip("{}"), dtype="float", sep=",")
+
+    parsed_tags = {}
+    for key, value in tags.items():
+        if value.startswith("{") and value.endswith("}"):
+            new_val = parsevec(value)
+            value = new_val if len(new_val) else value
+        else:
+            try:
+                value = int(value)
+            except (TypeError, ValueError):
+                try:
+                    value = float(value)
+                except (TypeError, ValueError):
+                    pass
+        parsed_tags[key] = value
+    return parsed_tags
+
+
 def open_rasterio(
     filename, parse_coordinates=None, chunks=None, cache=None, lock=None, masked=False
 ):
@@ -229,6 +250,8 @@ def open_rasterio(
     data : DataArray
         The newly created DataArray.
     """
+    parse_coordinates = True if parse_coordinates is None else parse_coordinates
+
     import rasterio
     from rasterio.vrt import WarpedVRT
 
@@ -255,6 +278,22 @@ def open_rasterio(
 
     manager = CachingFileManager(rasterio.open, filename, lock=lock, mode="r")
     riods = manager.acquire()
+
+    # open the subdatasets if they exist
+    if riods.subdatasets:
+        data_arrays = {}
+        for iii, subdataset in enumerate(riods.subdatasets):
+            rioda = open_rasterio(
+                subdataset,
+                parse_coordinates=iii == 0 and parse_coordinates,
+                chunks=chunks,
+                cache=cache,
+                lock=lock,
+                masked=masked,
+            )
+            data_arrays[rioda.name] = rioda
+        return Dataset(data_arrays)
+
     if vrt_params is not None:
         riods = WarpedVRT(riods, **vrt_params)
 
@@ -274,11 +313,10 @@ def open_rasterio(
     else:
         transform = riods.transform
 
-    parse = True if (parse_coordinates is None) else parse_coordinates
-    if transform.is_rectilinear and parse:
+    if transform.is_rectilinear and parse_coordinates:
         # 1d coordinates
         coords.update(affine_to_coords(riods.transform, riods.width, riods.height))
-    elif parse:
+    elif parse_coordinates:
         # 2d coordinates
         warnings.warn(
             "The file coordinates' transformation isn't "
@@ -290,7 +328,7 @@ def open_rasterio(
         )
 
     # Attributes
-    attrs = dict()
+    attrs = _parse_tags(riods.tags(1))
     encoding = dict()
     # Affine transformation matrix (always available)
     # This describes coefficients mapping pixel coordinates to CRS
@@ -341,7 +379,10 @@ def open_rasterio(
     if cache and chunks is None:
         data = indexing.MemoryCachedArray(data)
 
-    result = DataArray(data=data, dims=("band", "y", "x"), coords=coords, attrs=attrs)
+    da_name = attrs.pop("NETCDF_VARNAME", None)
+    result = DataArray(
+        data=data, dims=("band", "y", "x"), coords=coords, attrs=attrs, name=da_name
+    )
     result.encoding = encoding
 
     if hasattr(riods, "crs") and riods.crs:
diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py
index 8f30b619..077148da 100644
--- a/rioxarray/rioxarray.py
+++ b/rioxarray/rioxarray.py
@@ -91,7 +91,7 @@ def _generate_attrs(src_data_array, dst_affine, dst_nodata):
         new_attrs["_FillValue"] = fill_value
 
     # add raster spatial information
-    new_attrs["transform"] = tuple(dst_affine)
+    new_attrs["transform"] = tuple(dst_affine)[:6]
     new_attrs["grid_mapping"] = _get_grid_map_name(src_data_array)
 
     return new_attrs
diff --git a/sphinx/history.rst b/sphinx/history.rst
index 73e9c84b..b3daa8dc 100644
--- a/sphinx/history.rst
+++ b/sphinx/history.rst
@@ -1,6 +1,10 @@
 History
 =======
 
+0.0.10
+-----
+- Add support for opening netcdf/hdf files with `rioxarray.open_rasterio` (issue #32)
+
 0.0.9
 -----
 - Add `rioxarray.open_rasterio` (issue #7)


Test patch:
diff --git a/test/conftest.py b/test/conftest.py
index 8522e8db..63ca75aa 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -17,15 +17,17 @@ def _assert_xarrays_equal(input_xarray, compare_xarray, precision=7):
             assert "creation_date" in input_xr.attrs
 
         for attr in compare_xr.attrs:
-            if (
+            if attr == "transform":
+                assert_almost_equal(input_xr.attrs[attr], compare_xr.attrs[attr][:6])
+            elif (
                 attr != "_FillValue"
                 and attr not in UNWANTED_RIO_ATTRS
                 and attr != "creation_date"
             ):
                 try:
-                    assert input_xr.attrs[attr] == compare_xr.attrs[attr]
-                except ValueError:
                     assert_almost_equal(input_xr.attrs[attr], compare_xr.attrs[attr])
+                except (TypeError, ValueError):
+                    assert input_xr.attrs[attr] == compare_xr.attrs[attr]
 
     assert_attrs_equal(input_xarray, compare_xarray)
     if hasattr(input_xarray, "variables"):
diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py
index c9af0ba4..e0beae15 100644
--- a/test/integration/test_integration_rioxarray.py
+++ b/test/integration/test_integration_rioxarray.py
@@ -78,7 +78,9 @@ def interpolate_na_veris():
     )
 
 
-@pytest.fixture(params=[xarray.open_dataset, xarray.open_dataarray])
+@pytest.fixture(
+    params=[xarray.open_dataarray, rioxarray.open_rasterio, xarray.open_dataset]
+)
 def interpolate_na_nan(request):
     return dict(
         input=os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"),
@@ -87,7 +89,9 @@ def interpolate_na_nan(request):
     )
 
 
-@pytest.fixture(params=[xarray.open_dataset, xarray.open_dataarray])
+@pytest.fixture(
+    params=[xarray.open_dataset, xarray.open_dataarray, rioxarray.open_rasterio]
+)
 def modis_reproject_match(request):
     return dict(
         input=os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"),
@@ -103,7 +107,7 @@ def _mod_attr(input_xr, attr, val=None, remove=False):
             _mod_attr(input_xr[var], attr, val=val, remove=remove)
     else:
         if remove:
-            del input_xr.attrs[attr]
+            input_xr.attrs.pop(attr, None)
         else:
             input_xr.attrs[attr] = val
 
@@ -118,7 +122,9 @@ def _del_attr(input_xr, attr):
     _mod_attr(input_xr, attr, remove=True)
 
 
-@pytest.fixture(params=[xarray.open_dataset, xarray.open_dataarray])
+@pytest.fixture(
+    params=[xarray.open_dataset, xarray.open_dataarray, rioxarray.open_rasterio]
+)
 def modis_clip(request, tmpdir):
     return dict(
         input=os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"),
@@ -300,21 +306,23 @@ def test_transform_bounds():
 
 
 def test_reproject(modis_reproject):
+    mask_args = (
+        dict(masked=False)
+        if modis_reproject["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=False)
+    )
     with modis_reproject["open"](
-        modis_reproject["input"], mask_and_scale=False
-    ) as mda, modis_reproject["open"](
-        modis_reproject["compare"], mask_and_scale=False
-    ) as mdc:
+        modis_reproject["input"], **mask_args
+    ) as mda, modis_reproject["open"](modis_reproject["compare"], **mask_args) as mdc:
         mds_repr = mda.rio.reproject(modis_reproject["to_proj"])
         # test
         _assert_xarrays_equal(mds_repr, mdc)
 
 
-def test_reproject_3d(modis_reproject_3d):
-    with xarray.open_dataset(
-        modis_reproject_3d["input"], mask_and_scale=False
-    ) as mda, xarray.open_dataset(
-        modis_reproject_3d["compare"], mask_and_scale=False
+@pytest.fixture(params=[xarray.open_rasterio, rioxarray.open_rasterio])
+def test_reproject_3d(request, modis_reproject_3d):
+    with request.param(modis_reproject_3d["input"]) as mda, request.param(
+        modis_reproject_3d["compare"]
     ) as mdc:
         mds_repr = mda.rio.reproject(modis_reproject_3d["to_proj"])
         # test
@@ -322,11 +330,14 @@ def test_reproject_3d(modis_reproject_3d):
 
 
 def test_reproject__grid_mapping(modis_reproject):
+    mask_args = (
+        dict(masked=False)
+        if modis_reproject["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=False)
+    )
     with modis_reproject["open"](
-        modis_reproject["input"], mask_and_scale=False
-    ) as mda, modis_reproject["open"](
-        modis_reproject["compare"], mask_and_scale=False
-    ) as mdc:
+        modis_reproject["input"], **mask_args
+    ) as mda, modis_reproject["open"](modis_reproject["compare"], **mask_args) as mdc:
 
         # remove 'crs' attribute and add grid mapping
         mda.coords["spatial_ref"] = 0
@@ -376,11 +387,14 @@ def test_reproject__no_transform(modis_reproject):
 
 
 def test_reproject__no_nodata(modis_reproject):
+    mask_args = (
+        dict(masked=False)
+        if modis_reproject["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=False)
+    )
     with modis_reproject["open"](
-        modis_reproject["input"], mask_and_scale=False
-    ) as mda, modis_reproject["open"](
-        modis_reproject["compare"], mask_and_scale=False
-    ) as mdc:
+        modis_reproject["input"], **mask_args
+    ) as mda, modis_reproject["open"](modis_reproject["compare"], **mask_args) as mdc:
         orig_fill = _get_attr(mda, "_FillValue")
         _del_attr(mda, "_FillValue")
         _del_attr(mda, "nodata")
@@ -424,12 +438,17 @@ def test_reproject__no_nodata_masked(modis_reproject):
 
 
 def test_reproject_match(modis_reproject_match):
+    mask_args = (
+        dict(masked=False)
+        if modis_reproject_match["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=False)
+    )
     with modis_reproject_match["open"](
-        modis_reproject_match["input"], mask_and_scale=False
+        modis_reproject_match["input"], **mask_args
     ) as mda, modis_reproject_match["open"](
-        modis_reproject_match["compare"], mask_and_scale=False
+        modis_reproject_match["compare"], **mask_args
     ) as mdc, xarray.open_dataarray(
-        modis_reproject_match["match"], mask_and_scale=False
+        modis_reproject_match["match"]
     ) as mdm:
         # reproject
         mds_repr = mda.rio.reproject_match(mdm)
@@ -438,10 +457,15 @@ def test_reproject_match(modis_reproject_match):
 
 
 def test_reproject_match__masked(modis_reproject_match):
+    mask_args = (
+        dict(masked=True)
+        if modis_reproject_match["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=True)
+    )
     with modis_reproject_match["open"](
-        modis_reproject_match["input"]
+        modis_reproject_match["input"], **mask_args
     ) as mda, modis_reproject_match["open"](
-        modis_reproject_match["compare"]
+        modis_reproject_match["compare"], **mask_args
     ) as mdc, xarray.open_dataarray(
         modis_reproject_match["match"]
     ) as mdm:
@@ -452,10 +476,15 @@ def test_reproject_match__masked(modis_reproject_match):
 
 
 def test_reproject_match__no_transform_nodata(modis_reproject_match):
+    mask_args = (
+        dict(masked=True)
+        if modis_reproject_match["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=True)
+    )
     with modis_reproject_match["open"](
-        modis_reproject_match["input"]
+        modis_reproject_match["input"], **mask_args
     ) as mda, modis_reproject_match["open"](
-        modis_reproject_match["compare"]
+        modis_reproject_match["compare"], **mask_args
     ) as mdc, xarray.open_dataarray(
         modis_reproject_match["match"]
     ) as mdm:
@@ -559,21 +588,22 @@ def test_make_coords__attr_trans(request, modis_reproject):
 
 
 def test_interpolate_na(interpolate_na):
+    mask_args = (
+        dict(masked=False)
+        if interpolate_na["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=False)
+    )
     with interpolate_na["open"](
-        interpolate_na["input"], mask_and_scale=False
-    ) as mda, interpolate_na["open"](
-        interpolate_na["compare"], mask_and_scale=False
-    ) as mdc:
+        interpolate_na["input"], **mask_args
+    ) as mda, interpolate_na["open"](interpolate_na["compare"], **mask_args) as mdc:
         interpolated_ds = mda.rio.interpolate_na()
         # test
         _assert_xarrays_equal(interpolated_ds, mdc)
 
 
 def test_interpolate_na_veris(interpolate_na_veris):
-    with xarray.open_dataset(
-        interpolate_na_veris["input"], mask_and_scale=False
-    ) as mda, xarray.open_dataset(
-        interpolate_na_veris["compare"], mask_and_scale=False
+    with xarray.open_dataset(interpolate_na_veris["input"]) as mda, xarray.open_dataset(
+        interpolate_na_veris["compare"]
     ) as mdc:
         interpolated_ds = mda.rio.interpolate_na()
         # test
@@ -581,10 +611,8 @@ def test_interpolate_na_veris(interpolate_na_veris):
 
 
 def test_interpolate_na_3d(interpolate_na_3d):
-    with xarray.open_dataset(
-        interpolate_na_3d["input"], mask_and_scale=False
-    ) as mda, xarray.open_dataset(
-        interpolate_na_3d["compare"], mask_and_scale=False
+    with xarray.open_dataset(interpolate_na_3d["input"]) as mda, xarray.open_dataset(
+        interpolate_na_3d["compare"]
     ) as mdc:
         interpolated_ds = mda.rio.interpolate_na()
         # test
@@ -592,10 +620,15 @@ def test_interpolate_na_3d(interpolate_na_3d):
 
 
 def test_interpolate_na__nodata_filled(interpolate_na_filled):
+    mask_args = (
+        dict(masked=False)
+        if interpolate_na_filled["open"].__name__ == "open_rasterio"
+        else dict(mask_and_scale=False)
+    )
     with interpolate_na_filled["open"](
-        interpolate_na_filled["input"], mask_and_scale=False
+        interpolate_na_filled["input"], **mask_args
     ) as mda, interpolate_na_filled["open"](
-        interpolate_na_filled["compare"], mask_and_scale=False
+        interpolate_na_filled["compare"], **mask_args
     ) as mdc:
         if hasattr(mda, "variables"):
             for var in mda.rio.vars:
@@ -609,9 +642,13 @@ def test_interpolate_na__nodata_filled(interpolate_na_filled):
 
 
 def test_interpolate_na__all_nodata(interpolate_na_nan):
+    rio_opened = interpolate_na_nan["open"].__name__ == "open_rasterio"
+    mask_args = dict(masked=True) if rio_opened else dict(mask_and_scale=True)
     with interpolate_na_nan["open"](
-        interpolate_na_nan["input"]
-    ) as mda, interpolate_na_nan["open"](interpolate_na_nan["compare"]) as mdc:
+        interpolate_na_nan["input"], **mask_args
+    ) as mda, interpolate_na_nan["open"](
+        interpolate_na_nan["compare"], **mask_args
+    ) as mdc:
         if hasattr(mda, "variables"):
             for var in mda.rio.vars:
                 mda[var].values[~numpy.isnan(mda[var].values)] = numpy.nan
@@ -619,6 +656,9 @@ def test_interpolate_na__all_nodata(interpolate_na_nan):
             mda.values[~numpy.isnan(mda.values)] = numpy.nan
 
         interpolated_ds = mda.rio.interpolate_na()
+        if rio_opened and "__xarray_dataarray_variable__" in mdc:
+            mdc = mdc["__xarray_dataarray_variable__"]
+            mdc.attrs.pop("coordinates")
         # test
         _assert_xarrays_equal(interpolated_ds, mdc)
 
@@ -636,12 +676,13 @@ def test_load_in_geographic_dimensions():
         assert mda.red.rio.crs.to_epsg() == 4326
 
 
-def test_geographic_reproject():
+@pytest.fixture(params=[xarray.open_dataset, rioxarray.open_rasterio])
+def test_geographic_reproject(request):
     sentinel_2_geographic = os.path.join(
         TEST_INPUT_DATA_DIR, "sentinel_2_L1C_geographic.nc"
     )
     sentinel_2_utm = os.path.join(TEST_COMPARE_DATA_DIR, "sentinel_2_L1C_utm.nc")
-    with xarray.open_dataset(sentinel_2_geographic) as mda, xarray.open_dataset(
+    with request.param(sentinel_2_geographic) as mda, request.param(
         sentinel_2_utm
     ) as mdc:
         mds_repr = mda.rio.reproject("+init=epsg:32721")
@@ -650,14 +691,15 @@ def test_geographic_reproject():
         _assert_xarrays_equal(mds_repr, mdc)
 
 
-def test_geographic_reproject__missing_nodata():
+@pytest.fixture(params=[xarray.open_dataset, rioxarray.open_rasterio])
+def test_geographic_reproject__missing_nodata(request):
     sentinel_2_geographic = os.path.join(
         TEST_INPUT_DATA_DIR, "sentinel_2_L1C_geographic.nc"
     )
     sentinel_2_utm = os.path.join(
         TEST_COMPARE_DATA_DIR, "sentinel_2_L1C_utm__auto_nodata.nc"
     )
-    with xarray.open_dataset(sentinel_2_geographic) as mda, xarray.open_dataset(
+    with request.param(sentinel_2_geographic) as mda, request.param(
         sentinel_2_utm
     ) as mdc:
         mda.red.attrs.pop("nodata")
@@ -668,14 +710,15 @@ def test_geographic_reproject__missing_nodata():
         _assert_xarrays_equal(mds_repr, mdc)
 
 
-def test_geographic_resample_integer():
+@pytest.fixture(params=[xarray.open_dataset, rioxarray.open_rasterio])
+def test_geographic_resample_integer(request):
     sentinel_2_geographic = os.path.join(
         TEST_INPUT_DATA_DIR, "sentinel_2_L1C_geographic.nc"
     )
     sentinel_2_interp = os.path.join(
         TEST_COMPARE_DATA_DIR, "sentinel_2_L1C_interpolate_na.nc"
     )
-    with xarray.open_dataset(sentinel_2_geographic) as mda, xarray.open_dataset(
+    with request.param(sentinel_2_geographic) as mda, request.param(
         sentinel_2_interp
     ) as mdc:
         mds_interp = mda.rio.interpolate_na()
@@ -684,12 +727,11 @@ def test_geographic_resample_integer():
         _assert_xarrays_equal(mds_interp, mdc)
 
 
-def test_to_raster(tmpdir):
+@pytest.fixture(params=[xarray.open_dataarray, rioxarray.open_rasterio])
+def test_to_raster(request, tmpdir):
     tmp_raster = tmpdir.join("modis_raster.tif")
     test_tags = {"test": "1"}
-    with xarray.open_dataarray(
-        os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc")
-    ) as mda:
+    with request.param(os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc")) as mda:
         mda.rio.to_raster(str(tmp_raster), tags=test_tags)
         xds = mda.copy()
 
@@ -703,11 +745,10 @@ def test_to_raster(tmpdir):
         assert rds.tags() == {"AREA_OR_POINT": "Area", **test_tags}
 
 
-def test_to_raster_3d(tmpdir):
+@pytest.fixture(params=[xarray.open_dataset, rioxarray.open_rasterio])
+def test_to_raster_3d(request, tmpdir):
     tmp_raster = tmpdir.join("planet_3d_raster.tif")
-    with xarray.open_dataset(
-        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc")
-    ) as mda:
+    with request.param(os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc")) as mda:
         xds = mda.green.fillna(mda.green.rio.encoded_nodata)
         xds.rio._nodata = mda.green.rio.encoded_nodata
         xds.rio.to_raster(str(tmp_raster))