Raster Processing

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

19.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)
client

Client

Client-9c6bdf3a-7730-11ee-819f-6a0ce35f57ad

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: /user/halemohammad@clarku.edu/proxy/8787/status

Cluster Info

LocalCluster

2f5be50b

Dashboard: /user/halemohammad@clarku.edu/proxy/8787/status	Workers: 4
Total threads: 4	Total memory: 32.00 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-59bb8467-6757-4817-840d-67d3558e190c

Comm: tcp://127.0.0.1:46169	Workers: 4
Dashboard: /user/halemohammad@clarku.edu/proxy/8787/status	Total threads: 4
Started: Just now	Total memory: 32.00 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:41889	Total threads: 1
Dashboard: /user/halemohammad@clarku.edu/proxy/34163/status	Memory: 8.00 GiB
Nanny: tcp://127.0.0.1:33283
Local directory: /tmp/dask-scratch-space/worker-dsd87al7

Worker: 1

Comm: tcp://127.0.0.1:33807	Total threads: 1
Dashboard: /user/halemohammad@clarku.edu/proxy/34771/status	Memory: 8.00 GiB
Nanny: tcp://127.0.0.1:34591
Local directory: /tmp/dask-scratch-space/worker-yqd1tvmm

Worker: 2

Comm: tcp://127.0.0.1:45103	Total threads: 1
Dashboard: /user/halemohammad@clarku.edu/proxy/44105/status	Memory: 8.00 GiB
Nanny: tcp://127.0.0.1:35407
Local directory: /tmp/dask-scratch-space/worker-dfywurov

Worker: 3

Comm: tcp://127.0.0.1:38731	Total threads: 1
Dashboard: /user/halemohammad@clarku.edu/proxy/41753/status	Memory: 8.00 GiB
Nanny: tcp://127.0.0.1:35353
Local directory: /tmp/dask-scratch-space/worker-vhbpbw9f

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 2 -r 2
green_band = rioxarray.open_rasterio(sample_scene)
green_band.mean()

13.6 s ± 1.82 s per loop (mean ± std. dev. of 2 runs, 2 loops each)

Now, let’s add Dask.

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

14.3 s ± 3.93 s per loop (mean ± std. dev. of 2 runs, 2 loops each)

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

9.98 s ± 71.9 ms per loop (mean ± std. dev. of 2 runs, 2 loops each)

19.3. Mosaicking Rasters#

Let’s say you are interested to generate an RGB image of Sentinel-2 across your area of interest defined in the following bbox. We have defined our geomtry in three different formats suing bbox, aoi, and clipping_geom which will be used in various parts in the following. But they all refer to the same object.

from shapely.geometry import box
import shapely
import json

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

shapely.to_geojson(aoi_box)

'{"type":"Polygon","coordinates":[[[-69.46303157927986,41.01426426532743],[-69.46303157927986,42.0218426094371],[-70.74295109234716,42.0218426094371],[-70.74295109234716,41.01426426532743],[-69.46303157927986,41.01426426532743]]]}'

json.loads(shapely.to_geojson(aoi_box))

{'type': 'Polygon',
 'coordinates': [[[-69.46303157927986, 41.01426426532743],
   [-69.46303157927986, 42.0218426094371],
   [-70.74295109234716, 42.0218426094371],
   [-70.74295109234716, 41.01426426532743],
   [-69.46303157927986, 41.01426426532743]]]}

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

string indices must be integers, not 'str'

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

from pystac_client import Client

api_url = "https://earth-search.aws.element84.com/v1"
client = Client.open(api_url)
search_results = client.search(
    collections=["sentinel-2-l2a"],
    bbox = aoi_box.bounds,
    sortby=["-properties.datetime"],
    max_items=4,
)
items = search_results.item_collection()

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

string indices must be integers, not 'str'
string indices must be integers, not 'str'
string indices must be integers, not 'str'
string indices must be integers, not 'str'

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

We will use merge_datasets from rioxarray

import xarray as xr
import rioxarray
from rioxarray.merge import merge_datasets

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. 
    
    """
    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/DF/2023/10/S2B_19TDF_20231028_0_L2A/B02.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/DG/2023/10/S2B_19TDG_20231028_0_L2A/B02.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/CF/2023/10/S2A_19TCF_20231026_0_L2A/B02.tif',
 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/19/T/DF/2023/10/S2A_19TDF_20231026_0_L2A/B02.tif']

temp_dataarray = rioxarray.open_rasterio(red_urls[0])

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

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

red_datasets[0]

<xarray.Dataset>
Dimensions:      (band: 1, x: 10980, y: 10980)
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 4e+05 4e+05 4e+05 ... 5.097e+05 5.097e+05 5.098e+05
  * y            (y) float64 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06
    spatial_ref  int64 0
Data variables:
    red          (band, y, x) uint16 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 = [
    rioxarray.open_rasterio(url, lock=False, chunks=(1, 'auto', -1)).to_dataset(name="green") for url in green_urls
]

blue_datasets = [
    rioxarray.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)

Raster Processing

Contents

19. Raster Processing#

19.1. Reprojecting Raster#

19.2. Using `rioxarray` with Dask#

Client

Cluster Info

LocalCluster

Scheduler Info

Scheduler

Workers

Worker: 0

Worker: 1

Worker: 2

Worker: 3

19.3. Mosaicking Rasters#

Raster Processing

Contents

19. Raster Processing#

19.1. Reprojecting Raster#

19.2. Using rioxarray with Dask#

Client

Cluster Info

LocalCluster

Scheduler Info

Scheduler

Workers

Worker: 0

Worker: 1

Worker: 2

Worker: 3

19.3. Mosaicking Rasters#

19.2. Using `rioxarray` with Dask#