Raster Processing

15. Raster Processing#

The goal of this lecture is to learn how to do reprojections and mosaicking of rasters in Python. We will use rioxarray as the main package for these tasks.

15.2. Using `rioxarray` with Dask#

In this section, we will review how to use Dask with rioxarray to improve reading COG files. For maximum read performance, the chunking pattern you request with rio.open_rasterio should align with the internal chunking of the COG. Typically this means reading the data in a “row major” format: your chunks should be as wide as possible along the columns. (Check out this page for more information about this)

Let’s try this on a sample scene to read the raster, and calculate the mean value. We have a sample scene from green band of Sentinel-2 available on AWS S3:

import dask 
from dask.distributed import Client, LocalCluster
cluster = LocalCluster()
client = Client(cluster)
print(client.dashboard_link)

http://127.0.0.1:8787/status

sample_scene = "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/CF/2023/10/S2A_19TCF_20231026_0_L2A/B03.tif"

First, we use rioxarray without Dask:

%%timeit -n 1 -r 1
green_band = rxr.open_rasterio(sample_scene)
green_band.mean()

20.2 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)

Now, let’s add Dask.

%%timeit -n 1 -r 1
green_band_dask = rxr.open_rasterio(sample_scene, chunks=(1, 'auto', -1))
green_band_dask.mean().compute()

24.1 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)

%%timeit -n 1 -r 1
green_band_dask = rxr.open_rasterio(sample_scene, lock=False, chunks=(1, 'auto', -1))
green_band_dask.mean().compute()

12.6 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)

15.3. Mosaicking Rasters#

Let’s say you are interested in generating an RGB image of Sentinel-2 across your area of interest defined in the following aoi_box.

from shapely.geometry import box
import shapely
import json

aoi_box = box(-70.74295109234716, 41.01426426532743, -69.46303157927986, 42.0218426094371)

aoi_box

../_images/1b4395c77e09f16022e94a486fdf8ecf249d712b925c630e04c887fc035b9cb7.svg

First, let’s visualize this bbox

import leafmap

m = leafmap.Map(center=[41.633484, -70.154371], zoom=8, height="500px")
m.add_geojson(json.loads(shapely.to_geojson(aoi_box)))
m

We use the PySTAC client to search for imagery across this area:

import pystac_client
import planetary_computer

def retrieve_stac_items_bbox(url, collection, bbox, max_items = 10, pc_flag = False, query = None): 
    """
    This function retrieves the latest STAC items from a STAC API
    
    Parameters
    ----------
    url : str
        The STAC API URL
    collection : str
        Collection ID to search in the STAC API
    bbox : dict
       A tuple or list with coordinates of the target area for API search in the format (min x, min y, max x, max y)
    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. 
    max_items: integer
        Maximum number of items to retrieve from the STAC API. Default is 1. 

    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],
        bbox = bbox,
        query = query,
        sortby = ["-properties.datetime"],
        max_items = max_items
    )
    items = search_results.item_collection()

    return items

items = retrieve_stac_items_bbox(url = "https://earth-search.aws.element84.com/v1",
                                collection = "sentinel-2-l2a",
                                bbox = aoi_box.bounds,
                                max_items = 4,
                                pc_flag = False,
                                query = ["s2:nodata_pixel_percentage<1"]
                                )

for item in items:
    m.add_geojson(item.geometry, style = {"color": "red"})
m

So, as we can see in the map there are four scenes that overlap with our area of interest. We need to mosaick the overlapping Sentinel-2 scenes to generate one raster.

We will use merge_datasets from rioxarray

def gen_stac_asset_urls(items, asset):
    """
    This function receives an items collection returned by a STAC API, and returns
    the urls of the requested `asset` in a list. 

    Inputs:
        items : json collection
            A STAC items collection returned by STAC API
        asset : string
            Name of an asset present in the `items` collection

    Returns:
        urls : list
            List of all usls related to the `asset`
            
    """
    
    urls = []
    for item in items:
        urls.append(item.assets[asset].href)
    
    return urls

red_urls = gen_stac_asset_urls(items, "red")
green_urls = gen_stac_asset_urls(items, "green")
blue_urls = gen_stac_asset_urls(items, "blue")

blue_urls

['https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/CF/2024/10/S2A_19TCF_20241030_0_L2A/B02.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/DF/2024/10/S2A_19TDF_20241030_0_L2A/B02.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/CG/2024/10/S2A_19TCG_20241030_0_L2A/B02.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/DG/2024/10/S2A_19TDG_20241030_0_L2A/B02.tif']

temp_dataarray = rxr.open_rasterio(red_urls[0])

temp_dataset = rxr.open_rasterio(red_urls[0]).to_dataset(name="red")

red_datasets = [
    rxr.open_rasterio(url, lock=False, chunks=(1, 'auto', -1)).to_dataset(name="red") for url in red_urls
]

red_datasets[0]

<xarray.Dataset> Size: 241MB
Dimensions:      (band: 1, x: 10980, y: 10980)
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 88kB 3e+05 3e+05 3e+05 ... 4.098e+05 4.098e+05
  * y            (y) float64 88kB 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06
    spatial_ref  int64 8B 0
Data variables:
    red          (band, y, x) uint16 241MB dask.array<chunksize=(1, 6111, 10980), meta=np.ndarray>

mosaic_red = merge_datasets(red_datasets)

mosaic_red.red.rio.to_raster("red_mosaic.tif")

mosaic_red = mosaic_red.rio.clip([aoi_box], crs=4326)

mosaic_red.red.rio.to_raster("red_clipped.tif")

green_datasets = [
    rxr.open_rasterio(url, lock=False, chunks=(1, 'auto', -1)).to_dataset(name="green") for url in green_urls
]

blue_datasets = [
    rxr.open_rasterio(url, lock=False, chunks=(1, 'auto', -1)).to_dataset(name="blue") for url in blue_urls
]

mosaic_blue = merge_datasets(blue_datasets)
mosaic_blue = mosaic_blue.rio.clip([aoi_box], crs=4326)

mosaic_green = merge_datasets(green_datasets)
mosaic_green = mosaic_green.rio.clip([aoi_box], crs=4326)

rgb = xr.merge([mosaic_red, mosaic_green, mosaic_blue])

Raster Processing

Contents

15. Raster Processing#

15.1. Reprojecting a Raster File#

15.2. Using `rioxarray` with Dask#

15.3. Mosaicking Rasters#

Raster Processing

Contents

15. Raster Processing#

15.1. Reprojecting a Raster File#

15.2. Using rioxarray with Dask#

15.3. Mosaicking Rasters#

15.2. Using `rioxarray` with Dask#