Geospatial Raster Data in Python

18. 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.

18.1. Find a Raster 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:

from pystac_client import Client

api_url = "https://earth-search.aws.element84.com/v1"
client = Client.open(api_url)
collections = client.get_collections()
point = dict(type="Point", coordinates=(-71.822833, 42.250809))
search_results = client.search(
    collections=["sentinel-2-l2a"],
    intersects=point,
    query=["eo:cloud_cover<5", "s2:nodata_pixel_percentage>20"],
    sortby=["-properties.datetime"],
    max_items=5,
)
items = search_results.item_collection()

(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)

18.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.

18.2.1. Read the Metadata#

import rasterio
import os
import numpy as np

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

(Reminder: you can inspect the metadata of each item, and see what bands are available in the catalog using the assets.keys() property. e.g. items[0].assets.keys())

worcester_nir = rasterio.open(items[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 is a transform, check this reference.

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

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

All metada 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)}

18.2.2. 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, ..., 2784, 2754, 2574],
       [   0,    0,    0, ..., 2958, 2906, 2512],
       [   0,    0,    0, ..., 3000, 2814, 2354],
       ...,
       [   0,    0,    0, ...,  767,  758,  754],
       [   0,    0,    0, ...,  772,  738,  778],
       [   0,    0,    0, ...,  784,  794,  771]], dtype=uint16)

print(nir.dtype)

uint16

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

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

18.2.3. Visualize Band Data#

from rasterio import plot

plot.show(nir)

../_images/332031d1484e9fe7f5ca67b921460bf983ac5a85d91b2a47ca4d3348b783077b.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/c7eafc5e578337e2744843cbf8a702b31c7390441a17e46a35fd6c46ccc19872.png

<Axes: >

18.6. Writing Data#

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")

Task: Load the clipped raster and the original overview raster into QGIS and inspect their differences.

The difference is in their spatial resolution. Our clipped raster is for overview_level=1 which has a 40m spatial resolution. But the original oveerview has a native spatial resolution of 10m.

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

They represent the center coordinates of each pixel. You can verify this by comparing the bounds of the raster to the x and y coordinates of the top left pixel.

Geospatial Raster Data in Python

Contents

18. Geospatial Raster Data in Python#

18.1. Find a Raster Scene#

18.2. Working with Rasterio#

18.2.1. Read the Metadata#

18.2.2. Get Band Data#

18.2.3. Visualize Band Data#

18.3. Rioxarray#

18.4. Working with Multi-Band Data#

18.4.1. Clipping Raster#

18.5. Handling Large Rasters#

18.6. Writing Data#