How to use the rasterio.open function in rasterio

"""
    Create Cloud Optimized Geotiff.
    Parameters
    ----------
    src_path : str or PathLike object
        A dataset path or URL. Will be opened in "r" mode.
    This script is the rasterio equivalent of
    https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/validate_cloud_optimized_geotiff.py
    """
    errors = []
    details = {}

    if not GDALVersion.runtime().at_least("2.2"):
        raise Exception("GDAL 2.2 or above required")

    with rasterio.open(src_path) as src:
        if not src.driver == "GTiff":
            raise Exception("The file is not a GeoTIFF")

        filelist = [os.path.basename(f) for f in src.files]
        src_bname = os.path.basename(src_path)
        if len(filelist) > 1 and src_bname + ".ovr" in filelist:
            errors.append(
                "Overviews found in external .ovr file. They should be internal"
            )

        if src.width >= 512 or src.height >= 512:
            if not src.is_tiled:
                errors.append(
                    "The file is greater than 512xH or 512xW, but is not tiled"
                )

def get_projection(realpath, path):
    with rasterio.open(os.path.join(str(realpath), str(path))) as img:
        left, bottom, right, top = img.bounds
        spatial_reference = str(str(getattr(img, 'crs_wkt', None) or img.crs.wkt))
        geo_ref_points = {
            'ul': {'x': left, 'y': top},
            'ur': {'x': right, 'y': top},
            'll': {'x': left, 'y': bottom},
            'lr': {'x': right, 'y': bottom},
        }
        return geo_ref_points, spatial_reference

def from_geotiff(path_to_tif):
    with rasterio.drivers(CPL_DEBUG=True):
        with rasterio.open(path_to_tif) as src:
            b, g, r = src.read()

        total = np.zeros(r.shape, dtype=rasterio.uint16)
        for band in r, g, b:
            total += band
        total /= 3

    return total, src

import pprint

import rasterio
from rasterio.features import shapes

with rasterio.open('tests/data/shade.tif') as src:
    image = src.read(1)

# Print the first two shapes...
pprint.pprint(
    list(shapes(image))[:2]
)

def read_oneband_image_to_1dArray(image_path,nodata=None, ignore_small=None):

    if os.path.isfile(image_path) is False:
        raise IOError("error, file not exist: " + image_path)

    with rasterio.open(image_path) as img_obj:
        # read the all bands (only have one band)
        indexes = img_obj.indexes
        if len(indexes) != 1:
            raise IOError('error, only support one band')

        data = img_obj.read(indexes)
        data_1d = data.flatten()  # convert to one 1d, row first.

        if nodata is not None:
            data_1d = data_1d[data_1d != nodata]
        if ignore_small is not None:
            data_1d = data_1d[data_1d >= ignore_small ]

        return data_1d

def get_projection(path):
    with rasterio.open(str(path)) as img:
        left, bottom, right, top = img.bounds
        return {
            'spatial_reference': str(str(getattr(img, 'crs_wkt', None) or img.crs.wkt)),
            'geo_ref_points': {
                'ul': {
                    'x': left,
                    'y': top
                },
                'ur': {
                    'x': right,
                    'y': top
                },
                'll': {
                    'x': left,
                    'y': bottom
                },

valid_samples: array-like
        Numpy array of extracted raster values, typically 2d
    valid_coordinates: 2d array-like
        2D numpy array of xy coordinates of extracted values
    """

    # set the seed
    np.random.seed(seed=random_state)

    # determine number of GDAL rasters to sample
    n_features = len(rasters)

    # open rasters
    dataset = [0] * n_features
    for n in range(n_features):
        dataset[n] = rasterio.open(rasters[n], mode='r')

    # create np array to store randomly sampled data
    # we are starting with zero initial rows because data will be appended,
    # and number of columns are equal to n_features
    valid_samples = np.zeros((0, n_features))
    valid_coordinates = np.zeros((0, 2))

    # loop until target number of samples is satified
    satisfied = False

    n = size
    while satisfied is False:

        # generate random row and column indices
        Xsample = np.random.choice(range(0, dataset[0].shape[1]), n)
        Ysample = np.random.choice(range(0, dataset[0].shape[0]), n)

def estimateModel():
    
    # Open image files of predicted data and test data
    with rasterio.open(predicted_image_output_path, 'r') as prediction_dataset:      
        with rasterio.open(test_image_path, 'r') as test_labels_dataset:           
            
            #Find out the overlappin area of two images.
            #Because of tiling the prediction image is slightly smaller than the original clip.
            left = max(prediction_dataset.bounds.left,test_labels_dataset.bounds.left)
            bottom = max(prediction_dataset.bounds.bottom,test_labels_dataset.bounds.bottom)
            right = min(prediction_dataset.bounds.right,test_labels_dataset.bounds.right)
            top = min(prediction_dataset.bounds.top,test_labels_dataset.bounds.top)
            
            common_bbox = {
                        "type": "Polygon",
                        "coordinates": [[
                            [left, bottom],
                            [left, top],
                            [right, top],
                            [right, bottom],

crs = rasterio.crs.CRS.from_epsg(4326)  # WGS 84,
        x_left_top = xll_lon
        y_left_top = yll_lon + res*height
        profile = {'driver': 'GTiff', 'nodata': no_data, 'width': width, 'height': height, 'count': 1,
                   'crs': crs,
                   'transform': (x_left_top, res, 0.0, y_left_top , 0.0, -res),
                   }

        # And then change the band count to 1, set the
        # dtype to uint8, and specify LZW compression.
        profile.update(
            dtype=rasterio.float32,
            count=1,
            compress='lzw')

        with rasterio.open(output, 'w', **profile) as dst:
            dst.write(array_2d.astype(rasterio.float32), 1)

def get_model_coords(self):
        """Get LFP model coordinates for 1D and 2D mesh via BMI. The LFP model
        should be initialized first in order to access the variables."""

        logger.info('Getting LFP model coordinates.')

        i_ind, j_ind = np.where(np.logical_and(self.get_var('SGCwidth') > 0., self.get_var('DEM') != -9999))
        print i_ind.shape, j_ind.shape
        
        fn_map = join(self.model_data_dir,
                      self.model_config['dummysection']['DEMfile'])
        if not isfile(fn_map):
            raise IOError('landmask file not found')
        self._fn_landmask = fn_map

        with rasterio.open(fn_map, 'r') as ds:
            self.grid_index = ds.index
            self.model_grid_res = ds.res
            self.model_grid_bounds = ds.bounds
            self.model_grid_shape = ds.shape
            self.model_grid_transform = ds.transform

        list_x_coords, list_y_coords = self.model_grid_transform * (i_ind, j_ind)

        self.model_1d_coords = zip(list_x_coords, list_y_coords)
        self.model_1d_indices = zip(i_ind, j_ind)

        pass