Working with GeoTIFFs

Working with GeoTIFFs#

Prerequisites:

  • Required: None

  • Recommended: None

Introduction#

This is the first part of a series of tutorials to get users familiar with reVX.

reVX (NREL/reVX) is a tool used to support the reV model. It provides an interface to get data ready for reV modeling, as well as extracting and formatting reV outputs.

Most use cases of reVX would involve formatting geospatial raster and vector data for use within reV, as well as formatting reV outputs for uses outside reV. This tutorial will guide you through the process of working with GeoTIFF files using the Geotiff handler from the reVX library. GeoTIFF is a popular format for geospatial data, and the Geotiff handler provides a simple interface to load and extract various information from these files.

We’ll cover the following steps:

  1. Loading a GeoTIFF file using a context manager

  2. Retrieving the profile (metadata)

  3. Extracting band information

  4. Extracting data information

  5. Accessing the data values

  6. Extracting latitude and longitude coordinates

  7. Writing data to a GeoTIFF file

Let’s get started!

Downloading the data#

Before we dive into the code, we first have to download a sample TIFF from Siting Lab. In particular, we will be using data from [GDS24b].

If you have already downloaded the data, you can skip this step (just make sure path variables below are set correctly). We’ll start by defining the local file path destination:

DATA_FP = "airport_heliport_setbacks.tif"

Next, we can use a siting lab utility function to download the data. This function uses urllib (which is part of the Python standard library) under the hood.

Tip: You can view the contents of download_tiff_file online or by using the %load magic method in Jupyter. For the latter, simply type %load sl_utils at the top of an empty cell and execute it.
TIFF_URL = "https://data.openei.org/files/6120/airport_heliport_setbacks.tif"
download_tiff_file(TIFF_URL, DATA_FP, crop=True)
Note: The source TIFF files are large (90m resolution for all of CONUS), so we specified crop=True to crop the data immediately after downloading it to make it easier to work with. If you have a machine with sufficiently large memory (32GB+), or you are downloading the file in order to use it for analysis purposes, you should set crop=False.

Reading the data#

Now let’s use the Geotiff handler from the reVX library to open the TIFF file in the path defined.

Loading a GeoTIFF file using a context manager#

The following example will show the use of the handler within a context manager:

from reVX.handlers.geotiff import Geotiff

# We can use the Geotiff handler within a context manager
with Geotiff(DATA_FP) as geo:
    methods = [m for m in dir(geo) if not m.startswith("_")]

print("Geotiff methods:\n -", "\n - ".join(sorted(methods)))

Geotiff methods:
 - bands
 - close
 - dtype
 - iarr
 - lat_lon
 - latitude
 - longitude
 - meta
 - n_cols
 - n_rows
 - profile
 - shape
 - tiff_shape
 - values
 - write

Now we will use the Geotiff handler to inspect some properties of the file.

Retrieving the profile (metadata)#

We can use the profile attribute of the Geotiff class to get information on the profile:

with Geotiff(DATA_FP) as geo:
    profile = geo.profile
    print("GeoTIFF Profile:")
    pprint(profile)

GeoTIFF Profile:
{'blockxsize': 256,
 'blockysize': 256,
 'compress': 'lzma',
 'count': 1,
 'crs': '+init=epsg:5070',
 'driver': 'GTiff',
 'dtype': 'uint8',
 'height': 2000,
 'interleave': 'band',
 'nodata': 255.0,
 'tiled': True,
 'transform': (90.0, 0.0, 1829980.2632930684, 0.0, -90.0, 2297068.2309463923),
 'width': 2000}

Extracting band information#

We can extract the number of bands in the TIFF using the use the bands attribute:

with Geotiff(DATA_FP) as geo:
    bands = geo.bands
    print("Number of Bands:", bands)

Number of Bands: 1

Extracting data information#

We can extract information on the data in the TIFF using the following attributes:

with Geotiff(DATA_FP) as geo:
    # Determining the Data Type
    dtype = geo.dtype
    print("Data Type:", dtype)

    # Data shape
    shape = geo.shape
    print("Image shape:", shape)

Data Type: uint8
Image shape: (2000, 2000)

Accessing the data values#

We can extract the actual data from the TIFF using the values attribute:

with Geotiff(DATA_FP) as geo:
    # Extract data as a numpy array
    data_array = geo.values
    print(f"{type(data_array)} has a shape of {data_array.shape}")

<class 'numpy.ndarray'> has a shape of (1, 2000, 2000)

Extracting latitude and longitude coordinates#

So far we have used the Geotiff handler exclusively as a context manager. However, you may also use the Geotiff object without a context manager:

# Initialize a Geotiff object
geo = Geotiff(DATA_FP)

We can use the initialized Geotiff object to extract the (lat, lon) coordinates of all pixels using the lat_lon property:

# The lat_lon property returns the latitude and longitude values as a tuple
lat, lon = geo.lat_lon

print(f"latitude min, max: {lat.min(), lat.max()}")
print(f"longitude min, max: {lon.min(), lon.max()}")

latitude min, max: (39.7832, 41.729797)
longitude min, max: (-74.13154, -71.52781)

Alternatively, we can down-select to only the first 100x100 pixels (i.e. the Northwest corner of the raster) by requesting the "lat_lon" values from the handler:

# Like `geo.lat_lon`, this returns the latitude and longitude values as a tuple
lat, lon = geo["lat_lon", :100, :100]

print(f"NW Corner latitude min, max: {lat.min(), lat.max()}")
print(f"NW Corner longitude min, max: {lon.min(), lon.max()}")

NW Corner latitude min, max: (41.63378, 41.729797)
NW Corner longitude min, max: (-73.659775, -73.52976)

The latitude and longitude values can also be extracted individually by using the .latitude and `.longitude`` attributes of the handler, respectively.

lat = geo.latitude
lon = geo.longitude

print(f"latitude min, max: {lat.min(), lat.max()}")
print(f"longitude min, max: {lon.min(), lon.max()}")

latitude min, max: (39.7832, 41.729797)
longitude min, max: (-74.13154, -71.52781)

If using the handler without a context manager, use the close() method, to close the source object

geo.close()

Writing data#

We can also use the Geotiff handler to write data to a TIFF file. All we need to do is provide the data along with profile information, like so:

NEW_FILE = "gradient.tif"

profile = {
    "blockxsize": 256,
    "blockysize": 256,
    "compress": "lzma",
    "count": 1,
    "crs": "epsg:5070",
    "driver": "GTiff",
    "dtype": "uint8",
    "interleave": "band",
    "nodata": 255.0,
    "height": 2000,
    "width": 2000,
    "tiled": True,
    "transform": (90.0, 0.0, 1829980.263293, 0.0, -90.0, 2297068.230946),
}

new_data = np.arange(4000000).reshape(2000, 2000)

Geotiff.write(NEW_FILE, profile, new_data)

We can check that this operation succeeded:

with Geotiff(NEW_FILE) as geo:
    assert np.allclose(geo.values[0], new_data)

Conclusion#

In this tutorial, we have walked through the basic steps to load and explore GeoTIFF files using the Geotiff handler from the reVX library within a context manager. You should now be able to:

  • Retrieve metadata from a GeoTIFF file

  • Extract the values as a numpy array

  • Extract geographic coordinates (latitude and longitude)

  • Write data to a GeoTIFF file

Contents