How to use rioxarray - 10 common examples

def test_nc_attr_loading():
    with rioxarray.open_rasterio(
        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc")
    ) as rds:
        assert rds.dims == {"y": 10, "x": 10, "time": 2}
        assert rds.attrs == {"coordinates": "spatial_ref"}
        assert rds.y.attrs["units"] == "metre"
        assert rds.x.attrs["units"] == "metre"
        assert rds.time.encoding == {
            "units": "seconds since 2016-12-19T10:27:29.687763",
            "calendar": "proleptic_gregorian",
        }
        assert str(rds.time.values[0]) == "2016-12-19 10:27:29.687763"
        assert str(rds.time.values[1]) == "2016-12-29 12:52:42.347451"

def test_to_raster__dataset__mask_and_scale(tmpdir):
    output_raster = tmpdir.join("tmmx_20190121.tif")
    with rioxarray.open_rasterio(
        os.path.join(TEST_INPUT_DATA_DIR, "tmmx_20190121.nc")
    ) as rds:
        rds.isel(band=0).rio.to_raster(str(output_raster))

    with rioxarray.open_rasterio(str(output_raster)) as rdscompare:
        assert rdscompare.scale_factor == 0.1
        assert rdscompare.add_offset == 220.0
        assert rdscompare.long_name == "air_temperature"
        assert rdscompare.rio.crs == rds.rio.crs
        assert rdscompare.rio.nodata == rds.air_temperature.rio.nodata

def test_write_transform__from_read(tmp_path):
    xds = rioxarray.open_rasterio(
        os.path.join(TEST_COMPARE_DATA_DIR, "small_dem_3m_merged.tif"),
        parse_coordinates=False,
    )
    out_file = tmp_path / "test_geotransform.nc"
    xds.to_netcdf(out_file)
    xds2 = rioxarray.open_rasterio(out_file, parse_coordinates=False)
    assert_almost_equal(tuple(xds2.rio.transform()), tuple(xds.rio.transform()))
    assert xds.spatial_ref.GeoTransform == xds2.spatial_ref.GeoTransform

def test_missing_transform_bounds():
    xds = rioxarray.open_rasterio(
        os.path.join(TEST_COMPARE_DATA_DIR, "small_dem_3m_merged.tif"),
        parse_coordinates=False,
    )
    xds.coords["spatial_ref"].attrs.pop("GeoTransform")
    with pytest.raises(DimensionMissingCoordinateError):
        xds.rio.bounds()

def test_isel_window():
    with rioxarray.open_rasterio(
        os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc")
    ) as mda:
        assert (
            mda.rio.isel_window(Window.from_slices(slice(9, 12), slice(10, 12)))
            == mda.isel(x=slice(10, 12), y=slice(9, 12))
        ).all()

def test_merge_arrays__res():
    dem_test = os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc")
    with open_rasterio(dem_test, masked=True) as rds:
        rds.attrs = {
            "_FillValue": rds.rio.nodata,
            "grid_mapping": "spatial_ref",
            "crs": rds.attrs["crs"],
        }
        arrays = [
            rds.isel(x=slice(100), y=slice(100)),
            rds.isel(x=slice(100, 200), y=slice(100, 200)),
            rds.isel(x=slice(100), y=slice(100, 200)),
            rds.isel(x=slice(100, 200), y=slice(100)),
        ]
        merged = merge_arrays(arrays, res=(300, 300))

    assert_almost_equal(
        merged.rio.bounds(),
        (-7274009.649486291, 5003608.61015275, -7227509.649486291, 5050108.61015275),

def test_rasterio_vrt_with_src_crs(self):
        # Test open_rasterio() support of WarpedVRT with specified src_crs
        # create geotiff with no CRS and specify it manually
        with create_tmp_geotiff(crs=None) as (tmp_file, expected):
            src_crs = rasterio.crs.CRS.from_epsg(32618)
            with rasterio.open(tmp_file) as src:
                assert src.crs is None
                with rasterio.vrt.WarpedVRT(src, src_crs=src_crs) as vrt:
                    with rioxarray.open_rasterio(vrt) as da:
                        assert da.rio.crs == src_crs

def test_to_raster__dataset(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:
        mda.isel(time=0).rio.to_raster(str(tmp_raster))

    with rioxarray.open_rasterio(str(tmp_raster)) as rdscompare:
        assert rdscompare.scale_factor == 1.0
        assert rdscompare.add_offset == 0.0
        assert rdscompare.long_name == ("blue", "green")
        assert rdscompare.rio.crs == mda.rio.crs
        assert numpy.isnan(rdscompare.rio.nodata)

        (partial(rioxarray.open_rasterio, masked=True), False),
    ],
)
def test_to_raster_3d(open_method, windowed, tmpdir):
    tmp_raster = tmpdir.join("planet_3d_raster.tif")
    with open_method(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), windowed=windowed)
        xds_attrs = {
            key: str(value)
            for key, value in xds.attrs.items()
            if key not in ("add_offset", "nodata", "scale_factor", "transform")
        }

    with rasterio.open(str(tmp_raster)) as rds:
        assert rds.crs == xds.rio.crs

def test_open_group_filter__missing():
    with rioxarray.open_rasterio(
        os.path.join(TEST_INPUT_DATA_DIR, "PLANET_SCOPE_3D.nc"),
        variable="blue",
        group=["non-existent"],
    ) as rds:
        assert list(rds.data_vars) == []