<jupyter_start><jupyter_text>Read GeoTIFF data through a BMI This notebook describes how to open and read data from GeoTIFF filesusing a [Basic Model Interface](https://bmi.readthedocs.io/) (BMI) through the `BmiGeoTiff` class included in the `bmi-geotiff` package.  Setup To ensure all dependencies are met, set up a conda environment using the environment file found in the root directory of this repository:```conda env create --file=environment.yml```Then install the `bmi-geotiff` package:```make install``` Import a set of libraries for later use:<jupyter_code>import numpy as np
from rasterio.plot import show<jupyter_output><empty_output><jupyter_text>Open a file Import the `BmiGeoTiff` class from the `bmi-geotiff` package:<jupyter_code>from bmi_geotiff import BmiGeoTiff<jupyter_output><empty_output><jupyter_text>Create an instance of this class.<jupyter_code>m = BmiGeoTiff()<jupyter_output><empty_output><jupyter_text>Calling `help` on the instance displays all the BMI methods that are available.<jupyter_code>help(m)<jupyter_output>Help on BmiGeoTiff in module bmi_geotiff.bmi object:

class BmiGeoTiff(bmipy.bmi.Bmi)
 |  BmiGeoTiff() -> None
 |  
 |  BMI-mediated access to data and metadata in a GeoTIFF file.
 |  
 |  Method resolution order:
 |      BmiGeoTiff
 |      bmipy.bmi.Bmi
 |      abc.ABC
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self) -> None
 |      Initialize self.  See help(type(self)) for accurate signature.
 |  
 |  finalize(self) -> None
 |      Perform tear-down tasks for the model.
 |      
 |      Perform all tasks that take place after exiting the model's time
 |      loop. This typically includes deallocating memory, closing files and
 |      printing reports.
 |  
 |  get_component_name(self) -> str
 |      Name of the component.
 |      
 |      Returns
 |      -------
 |      str
 |          The name of the component.
 |  
 |  get_current_time(self) -> float
 |      Current time of the model.
 |      
 |      Returns
 |      -------
 |      float
 |          [...]<jupyter_text>The first step in using a BMI is calling the `initialize` method.This method requires a configuration file that provides initial values for the `GeoTiff` library wrapped by the BMI.A sample configuration file is provided in the current directory.<jupyter_code>ls
cat config.yaml<jupyter_output>bmi-geotiff:
  filename: RGB.byte.tif<jupyter_text>In this case, the configuration file simply lists the path to a GeoTIFF file(here, the test file RGB.byte.tif from the [rasterio](https://rasterio.readthedocs.io/) project).Call `initialize` with the sample configuration file:<jupyter_code>m.initialize("config.yaml")<jupyter_output><empty_output><jupyter_text>The GeoTIFF file listed in the configuration file has now been opened,and the information is contains can be accessed through BMI methods.  Access data through the BMI Now that we've opened the GeoTIFF file, let's access the data and metadata it contains through the BMI.This will take a few steps.It may seem cumbersome at first, but there's payoff at the end.Start by displying the names of the variables exposed through the BMI.<jupyter_code>m.get_output_var_names()<jupyter_output><empty_output><jupyter_text>The (long) names used for these variables are instances of [CSDMS Standard Names](https://csdms.colorado.edu/wiki/CSDMS_Standard_Names).Standard Names are intended to be unambiguous; the tradeoff is that they tend to be long.Here, the first variable is for the raster data stored in the file. Find the data type of the raster.<jupyter_code>dtype = m.get_var_type("gis__raster_data")
dtype<jupyter_output><empty_output><jupyter_text>Within the BMI, functions that describe the grids that variables are defined on take an index instead of a variable name.Get the grid index for the raster data variable.<jupyter_code>grid = m.get_var_grid("gis__raster_data")
grid<jupyter_output><empty_output><jupyter_text>Then find the total size of the raster data.<jupyter_code>size = m.get_grid_size(grid)
size<jupyter_output><empty_output><jupyter_text>Next, get the raster data values.Two notes, however:* As a rule, memory should not be allocated within a BMI. This leads to the un-Pythonic way that we get the data--first creating an empty array, then passing it to a BMI function to receive values.* BMI arrays are flattened. This obviates array ordering issues between languages, but it does make >1D data harder to work with.Allocate an array for the raster data.<jupyter_code>raster = np.ndarray(size, dtype)
raster<jupyter_output><empty_output><jupyter_text>Get the data.<jupyter_code>m.get_value("gis__raster_data", raster)<jupyter_output><empty_output><jupyter_text>Note that the array is one-dimensional.<jupyter_code>raster.shape<jupyter_output><empty_output><jupyter_text>Reshape the data Like all BMI arrays, the raster values returned from the BMI `get_value` function are flattened.Let's restore their original dimensionality. First, determine the dimensionality of the raster variable.<jupyter_code>rank = m.get_grid_rank(grid)
rank<jupyter_output><empty_output><jupyter_text>Get the dimensions of the raster data, first creating an array to store their values.<jupyter_code>shape = np.empty(rank, dtype=int)
m.get_grid_shape(grid, shape)<jupyter_output><empty_output><jupyter_text>Reshape the raster data, creating a new array.<jupyter_code>rasterRGB = raster.reshape(shape)
rasterRGB.shape<jupyter_output><empty_output><jupyter_text>Get map projection information The data in the GeoTIFF file are georeferenced.The second and third variables exposed through the BMI,"gis__coordinate_reference_system" and "gis__gdal_geotransform", respectively,contain the proj.4 string and the [GDAL geotransform](https://gdal.org/user/raster_data_model.htmlaffine-geotransform) for the data.Get the geotransform through the BMI, following a process similar to what was used to obtain the raster data:* get variable type* get id of grid on which variable is defined* get size of grid* use the above to allocate an array for the transform* get the transform<jupyter_code>dtype = m.get_var_type("gis__gdal_geotransform")
dtype
grid = m.get_var_grid("gis__gdal_geotransform")
grid
size = m.get_grid_size(grid)
size
transform = np.ndarray(size, dtype)
m.get_value("gis__gdal_geotransform", transform)
transform<jupyter_output><empty_output><jupyter_text>Visualize Let's visualize the raster data as an image, with a little help from rasterio.<jupyter_code>show(rasterRGB, transform=list(transform))<jupyter_output><empty_output><jupyter_text>Conclusion Last, call the BMI `finalize` function.<jupyter_code>m.finalize()<jupyter_output><empty_output>