How to use the xarray.open_rasterio function in xarray

def test_utm(self):
        with create_tmp_geotiff() as (tmp_file, expected):
            with xr.open_rasterio(tmp_file) as rioda:
                assert_allclose(rioda, expected)
                assert rioda.attrs["scales"] == (1.0, 1.0, 1.0)
                assert rioda.attrs["offsets"] == (0.0, 0.0, 0.0)
                assert rioda.attrs["descriptions"] == ("d1", "d2", "d3")
                assert rioda.attrs["units"] == ("u1", "u2", "u3")
                assert isinstance(rioda.attrs["crs"], str)
                assert isinstance(rioda.attrs["res"], tuple)
                assert isinstance(rioda.attrs["is_tiled"], np.uint8)
                assert isinstance(rioda.rio._cached_transform(), Affine)
                np.testing.assert_array_equal(
                    rioda.attrs["nodatavals"], [np.NaN, np.NaN, np.NaN]
                )

            # Check no parse coords
            with xr.open_rasterio(tmp_file, parse_coordinates=False) as rioda:
                assert "x" not in rioda.coords

def test_rasterio_environment(self):
        with create_tmp_geotiff() as (tmp_file, expected):
            # Should fail with error since suffix not allowed
            with pytest.raises(Exception):
                with rasterio.Env(GDAL_SKIP="GTiff"):
                    with xr.open_rasterio(tmp_file) as actual:
                        assert_allclose(actual, expected)

def test_no_mftime(self):
        # rasterio can accept "filename" urguments that are actually urls,
        # including paths to remote files.
        # In issue #1816, we found that these caused dask to break, because
        # the modification time was used to determine the dask token. This
        # tests ensure we can still chunk such files when reading with
        # rasterio.
        with create_tmp_geotiff(
            8, 10, 3, transform_args=[1, 2, 0.5, 2.0], crs="+proj=latlong"
        ) as (tmp_file, expected):
            with mock.patch("os.path.getmtime", side_effect=OSError):
                with xr.open_rasterio(tmp_file, chunks=(1, 2, 2)) as actual:
                    assert isinstance(actual.data, da.Array)
                    assert_allclose(actual, expected)

def readImage(image_folder_fp):
    print("Reading Sentinel image from: %s" % (image_folder_fp))

    ### Rather than figuring out what the filepath inside SAFE folder is, this is just finding the red and nir files with correct endings
    for subdir, dirs, files in os.walk(image_folder_fp):
        for file in files:
            if file.endswith("_B04_10m.jp2"):
                red_fp = os.path.join(subdir,file)
            if file.endswith("_B08_10m.jp2"):
                nir_fp = os.path.join(subdir,file)

    ### Read the red and nir band files to xarray and with the chunk-option to dask
    red = xr.open_rasterio(red_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})
    nir = xr.open_rasterio(nir_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})

    ### Scale the image values back to real reflectance values
    red = red /10000
    nir = nir /10000

    return red,nir

def read(self):
        """Read the image."""
        dataset = rasterio.open(self.finfo['filename'])

        # Create area definition
        if hasattr(dataset, 'crs') and dataset.crs is not None:
            self.area = utils.get_area_def_from_raster(dataset)

        data = xr.open_rasterio(dataset, chunks=(1, CHUNK_SIZE, CHUNK_SIZE))
        attrs = data.attrs.copy()

        # Rename to Satpy convention
        data = data.rename({'band': 'bands'})

        # Rename bands to [R, G, B, A], or a subset of those
        data['bands'] = BANDS[data.bands.size]

        # Mask data if alpha channel is present
        try:
            data = mask_image_data(data)
        except ValueError as err:
            logger.warning(err)

        data.attrs = attrs
        self.file_content['image'] = data

S=6205900.
W=344520.
E=344620.
plot_bbox = np.asarray([[W,N],[E,N],[E,S],[W,S]])
pts, starting_ids, trees = io.load_lidar_data_by_polygon(las_list,plot_bbox,laz_files=False,max_pts_per_tree = 5*10**5)
N_trees = len(trees)
# filter LiDAR data by return classification
pts[np.any((pts[:,4]==3,pts[:,4]==4,pts[:,4]==5),axis=0),4]=1
pts[pts[:,2]<0,2]=0

# Load sensitivity analysis to spatial resolution
resolution_sensitivity = np.load('SEOS_MH_sensitivity_resolution.npy')[()]
density_sensitivity = np.load('SEOS_MH_sensitivity_pulse_density.npy')[()]

# Load 3D maps of canopy density
pad = xr.open_rasterio(pad_file)

# other parameters
max_height = 40.
layer_thickness = 0.5
heights = np.arange(0.,max_height,layer_thickness)+0.5
n_layers = heights.size
plot_width = 100.
sample_res = 5.
kappa = 1.
n_iter = 100

# bin lidar returns
heights,first_return_profile,n_ground_returns = LAD.bin_returns(pts, max_height, layer_thickness)

"""
PART B: Plot canopy profiles

def open(self, path, fmt, with_nodata=False, isel_band=None):
        "Open raster and return in requested format"

        if fmt is None or fmt.lower() == 'xarray':
            xarr = xr.open_rasterio(path, parse_coordinates=True)
            if isel_band is not None:
                xarr = xarr.isel(band=0)
            if with_nodata:
                xarr = convert_nodata_to_nans(xarr)
            return xarr

        raster_data = self.read(path)
        if fmt.lower() == 'rasterio':
            return raster_data

        if fmt.lower() == 'array':
            return raster_data.read()

        raise NotImplementedError('format %s not recognized' % fmt)

def readImage(image_folder_fp):
    print("Reading Sentinel image from: %s" % (image_folder_fp))

    ### Rather than figuring out what the filepath inside SAFE folder is, this is just finding the red and nir files with correct endings
    for subdir, dirs, files in os.walk(image_folder_fp):
        for file in files:
            if file.endswith("_B04_10m.jp2"):
                red_fp = os.path.join(subdir,file)
            if file.endswith("_B08_10m.jp2"):
                nir_fp = os.path.join(subdir,file)

    ### Read the red and nir band files to xarray and with the chunk-option to dask
    red = xr.open_rasterio(red_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})
    nir = xr.open_rasterio(nir_fp, chunks={'band': 1, 'x': 1024, 'y': 1024})

    ### Scale the image values back to real reflectance values
    red = red /10000
    nir = nir /10000

    return red,nir

Using rasterio's projection information for more accurate plots.

This example extends :ref:`recipes.rasterio` and plots the image in the
original map projection instead of relying on pcolormesh and a map
transformation.
"""

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

import xarray as xr

# Read the data
url = "https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif"
da = xr.open_rasterio(url)

# The data is in UTM projection. We have to set it manually until
# https://github.com/SciTools/cartopy/issues/813 is implemented
crs = ccrs.UTM("18N")

# Plot on a map
ax = plt.subplot(projection=crs)
da.plot.imshow(ax=ax, rgb="band", transform=crs)
ax.coastlines("10m", color="r")
plt.show()