Geospatial Raster Data in Python

14. Geospatial Raster Data in Python#

Attribution: Parts of this tutorial are developed based on the content from the following great sources: Parts of this tutorial are developed based on the content from the following great sources: Introduction to Python for Geographic Data AnalysisIntroduction to Python for Geographic Data Analysis; Introduction to Raster Data in PythonIntroduction to Raster Data in Python; and Introduction to Geospatial Raster and Vector Data with PythonIntroduction to Geospatial Raster and Vector Data with Python

The purpose of this lecture is to learn rasterio and rioxarray for loading and working with raster data.

You can access this notebook (in a Docker image) on this GitHub repo.

14.1. Find a Satellite Scene#

First, Let’s find a Sentinel-2 scene to work with. We will use the PySTAC client to retrieve the latest cloud free scene across Clark University campus:

import pystac_client
import planetary_computer

def retrieve_latest_stac_item(url, collection, point, pc_flag = False, query = None): 
    """
    This function retrieves the latest STAC item from a STAC API
    
    Parameters
    ----------
    url : str
        The STAC API URL
    collection : str
        Collection ID to search in the STAC API
    point : dict
       A GeoJOSN Point feature with coordinates of the target point for API search
    pc_flag : boolean
        A boolean flag to specify whether the STAC API is Microsoft's planetary computer API. Default is set to False.
    query : list
        A list of query parameters for STAC API. Default is set to None. 

    Returns
    -------
    items : Generator
        The STAC item that has the latest date of acquisition given all the query parameters. 
        
    """
    if pc_flag:
        modifier = planetary_computer.sign_inplace
    else:
        modifier = None
        
    catalog = pystac_client.Client.open(
        url = url,
        modifier = modifier,
    )

    search_results = catalog.search(
        collections = [collection],
        intersects = point,
        query = query,
        sortby = ["-properties.datetime"],
        max_items = 1
    )
    item = search_results.item_collection()

    return item

# A point in Clark campus:
point = dict(type="Point", coordinates=(-71.822833, 42.250809))

# Query parameters for low cloud cover and high missing pixels.
query = ["eo:cloud_cover<5", "s2:nodata_pixel_percentage>20"]

item = retrieve_latest_stac_item(url = "https://earth-search.aws.element84.com/v1",
                                 collection = "sentinel-2-l2a",
                                 point = point, 
                                 query = query
                                )      

(Note that we also added s2:nodata_pixel_percentage to our query parameters to find scenes that have some missing pixels. We want to use these to demonstrate how you would handle missing data later on)

14.2. Working with Rasterio#

Rasterio is a highly useful module for raster processing which you can use for reading and writing several different raster formats in Python. Rasterio is based on GDAL and Python automatically registers all known GDAL drivers for reading supported formats when importing the module.

14.2.1. Read the Metadata#

import rasterio
import os
import numpy as np

Let’s open the nir band of the first item:

worcester_nir = rasterio.open(item[0].assets["nir"].href)

And let’s see what is the type of worcester_nir:

type(worcester_nir)

rasterio.io.DatasetReader

You can see that our raster variable is a rasterio.io.DatasetReader type which means that we have opened the file for reading. This does not contain the value of the bands yet.

You can retrieve the metadata of this raster as following:

print("CRS of the scene is:", worcester_nir.crs)
print("Transform of the scene is:", worcester_nir.transform)
print("Width of the scene is:", worcester_nir.width)
print("Height of the scene is:", worcester_nir.height)
print("Number of bands in the scene is:", worcester_nir.count)
print("Bounds of the scene are:", worcester_nir.bounds)
print("The no data values for all channels are:", worcester_nir.nodatavals)

CRS of the scene is: EPSG:32618
Transform of the scene is: | 10.00, 0.00, 699960.00|
| 0.00,-10.00, 4700040.00|
| 0.00, 0.00, 1.00|
Width of the scene is: 10980
Height of the scene is: 10980
Number of bands in the scene is: 1
Bounds of the scene are: BoundingBox(left=699960.0, bottom=4590240.0, right=809760.0, top=4700040.0)
The no data values for all channels are: (0.0,)

If you need a reminder about what transform is, check this reference.

Alternatively, you can retrieve all metadata of the scene using the following:

print("All metadata of the scene:")
print(worcester_nir.meta)

All metadata of the scene:
{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': 0.0, 'width': 10980, 'height': 10980, 'count': 1, 'crs': CRS.from_epsg(32618), 'transform': Affine(10.0, 0.0, 699960.0,
       0.0, -10.0, 4700040.0)}

14.2.2. Exercise 1#

You can retrieve the same metadata from the STAC item properties. Write a code that prints the same metadata as printed in the previous cell but using the STAC item properties:

14.2.3. Get Band Data#

You can read the data from one or multiple bands using the read method:

nir = worcester_nir.read(1)

Note: rasterio band numbering starts from 1 unlike Python that starts from 0.

Let’s examine the data type this band:

nir

array([[   0,    0,    0, ..., 2750, 3100, 3396],
       [   0,    0,    0, ..., 3392, 4028, 3228],
       [   0,    0,    0, ..., 3608, 3828, 3132],
       ...,
       [   0,    0,    0, ...,    1,    1,    1],
       [   0,    0,    0, ...,   11,    6,    8],
       [   0,    0,    0, ...,    2,    1,    1]], dtype=uint16)

print(nir.dtype)

uint16

This data type is a 16-bit unsigned integer. Unsigned integer is always equal to or greater than zero but signed integer can store negative values too. For example, an unsigned 16-bit integer can store 2^16 (=65,536) values ranging from 0 to 65,535.

14.2.4. Visualize Band Data#

Let’s now visualize this band. We can use the plot function from rasterio:

from rasterio import plot

plot.show(nir)

../_images/49d6f4f246b264e2fe8a9489b302cc397a56ce7235fa2a8034161732dc23ca7a.png

<Axes: >

As you can see in this case the band data is treated as a regular (non-geospatial) array and the axes start from 0 instead of the actual coordinates of the raster. This happens as the nir variable is a numpy array and has no geospatial attirbutes. You can use the rasterio dataset worcester_nir and get the plot with the actual bounds:

plot.show(worcester_nir)

../_images/1b6eb29a5860ecb8c16751cca97d53224aaf3b65468792f07d79052cfd577554.png

<Axes: >

14.7. Writing Data to Disk#

You can use the .rio.to_raster function to write a rioxarray.DataArray to disk. Let’s write the clipped version of our multi-band imagery to disk:

worcester_clipped = worcester_overview.rio.clip([clipping_geom])

worcester_clipped.rio.to_raster("worcester_clipped.tif")

14.8. Exercise 2#

Question: What do x and y represent? Are they the center coordinates of each pixel, or the edge?

Geospatial Raster Data in Python

Contents

14. Geospatial Raster Data in Python#

14.1. Find a Satellite Scene#

14.2. Working with Rasterio#

14.2.1. Read the Metadata#

14.2.2. Exercise 1#

14.2.3. Get Band Data#

14.2.4. Visualize Band Data#

14.3. Working with Rioxarray#

14.4. Working with Multi-Band Data#

14.5. Clipping a Raster#

14.6. Handling Large Rasters#

14.7. Writing Data to Disk#

14.8. Exercise 2#