Vector segmentation masks#

Clouds float by, water flows on; in movement there is no grasping, in Chan there is no settling

For 🧑‍🏫 supervised machine learning, labels 🏷️ are needed in addition to the input image 🖼️. Here, we’ll step through an example workflow on matching vector 🚏 label data (points, lines, polygons) to 🛰️ Earth Observation data inputs. Specifically, this tutorial will cover:

  • Reading shapefiles 📁 directly from the web via pyogrio

  • Rasterizing vector polygons from a geopandas.GeoDataFrame to an xarray.DataArray

  • Pairing 🛰️ satellite images with the rasterized label masks and feeding them into a DataLoader

🎉 Getting started#

These are the tools 🛠️ you’ll need.

import matplotlib.pyplot as plt
import numpy as np
import planetary_computer
import pyogrio
import pystac
import torch
import torchdata
import xarray as xr
import zen3geo

0️⃣ Find cloud-hosted raster and vector data ⛳#

In this case study, we’ll look at the flood water extent over Johor, Malaysia 🇲🇾 on 15 Dec 2019 that were digitized by 🇺🇳 UNITAR-UNOSAT’s rapid mapping service over Synthetic Aperture Radar (SAR) 🛰️ images. Specifically, we’ll be using the 🇪🇺 Sentinel-1 Ground Range Detected (GRD) product’s VV polarization channel.

🔗 Links:

To start, let’s get the 🛰️ satellite scene we’ll be using for this tutorial.

item_url = "https://planetarycomputer.microsoft.com/api/stac/v1/collections/sentinel-1-grd/items/S1A_IW_GRDH_1SDV_20191215T224757_20191215T224822_030365_037955"

# Load the individual item metadata and sign the assets
item = pystac.Item.from_file(item_url)
signed_item = planetary_computer.sign(item)
signed_item

This is how the Sentinel-1 🩻 image looks like over Johor in Peninsular Malaysia on 15 Dec 2019.

Sentinel-1 image over Johor, Malaysia on 20191215

Load and reproject image data 🔄#

To keep things simple, we’ll load just the VV channel into a DataPipe via zen3geo.datapipes.RioXarrayReader (functional name: read_from_rioxarray) 😀.

url = signed_item.assets["vv"].href
dp = torchdata.datapipes.iter.IterableWrapper(iterable=[url])
# Reading lower resolution grid using overview_level=3
dp_rioxarray = dp.read_from_rioxarray(overview_level=3)
dp_rioxarray

RioXarrayReaderIterDataPipe

The Sentinel-1 image from Planetary Computer comes in longitude/latitude 🌐 geographic coordinates by default (OGC:CRS84). To make the pixels more equal 🔲 area, we can project it to a 🌏 local projected coordinate system instead.

def reproject_to_local_utm(dataarray: xr.DataArray, resolution: float=80.0) -> xr.DataArray:
    """
    Reproject an xarray.DataArray grid from OGC:CRS84 to a local UTM coordinate
    reference system.
    """
    # Estimate UTM coordinate reference from a single pixel
    pixel = dataarray.isel(y=slice(0, 1), x=slice(0,1))
    new_crs = dataarray.rio.reproject(dst_crs="OGC:CRS84").rio.estimate_utm_crs()

    return dataarray.rio.reproject(dst_crs=new_crs, resolution=resolution)

dp_reprojected = dp_rioxarray.map(fn=reproject_to_local_utm)

Note

Universal Transverse Mercator (UTM) isn’t actually an equal-area projection system. However, Sentinel-1 🛰️ satellite scenes from Copernicus are usually distributed in a UTM coordinate reference system, and UTM is typically a close enough 🤏 approximation to the local geographic area, or at least it won’t matter much when we’re looking at spatial resolutions over several 10s of metres 🙂.

Hint

For those wondering what OGC:CRS84 is, it is the longitude/latitude version of EPSG:4326 🌐 (latitude/longitude). I.e., it’s a matter of axis order, with OGC:CRS84 being x/y and EPSG:4326 being y/x.

🔖 References:

Transform and visualize raster data 🔎#

Let’s visualize 👀 the Sentinel-1 image, but before that, we’ll transform 🔄 the VV data from linear to decibel scale.

def linear_to_decibel(dataarray: xr.DataArray) -> xr.DataArray:
    """
    Transforming the input xarray.DataArray's VV or VH values from linear to
    decibel scale using the formula ``10 * log_10(x)``.
    """
    # Mask out areas with 0 so that np.log10 is not undefined
    da_linear = dataarray.where(cond=dataarray != 0)
    da_decibel = 10 * np.log10(da_linear)
    return da_decibel

dp_decibel = dp_reprojected.map(fn=linear_to_decibel)
dp_decibel

MapperIterDataPipe

As an aside, we’ll be using the Sentinel-1 image datapipe twice later, once as a template to create a blank canvas 🎞️, and another time by itself 🪞. This requires forking 🍴 the DataPipe into two branches, which can be achieved using torchdata.datapipes.iter.Forker (functional name: fork).

dp_decibel_canvas, dp_decibel_image = dp_decibel.fork(num_instances=2)
dp_decibel_canvas, dp_decibel_image

(_ChildDataPipe, _ChildDataPipe)

Now to visualize the transformed Sentinel-1 image 🖼️. Let’s zoom in 🔭 to one of the analysis extent areas we’ll be working on later.

it = iter(dp_decibel_image)
dataarray = next(it)

da_clip = dataarray.rio.clip_box(minx=371483, miny=190459, maxx=409684, maxy=229474)
da_clip.isel(band=0).plot.imshow(figsize=(11.5, 9), cmap="Blues_r", vmin=18, vmax=26)

<matplotlib.image.AxesImage at 0x7f91c72a2d90>
_images/7b303a6f550bf71b207256fd77cec4e60d00a5093bdfb87d61ef1d75ec391ce8.png

Notice how the darker blue areas 🔵 tend to correlate more with water features like the meandering rivers and the 🐚 sea on the NorthEast. This is because the SAR 🛰️ signal which is side looking reflects off flat water bodies like a mirror 🪞, with little energy getting reflected 🙅 back directly to the sensor (hence why it looks darker ⚫).

Load and visualize cloud-hosted vector files 💠#

Let’s now load some vector data from the web 🕸️. These are polygons of the segmented 🌊 water extent digitized by UNOSAT’s AI Based Rapid Mapping Service. We’ll be converting these vector polygons to 🌈 raster masks later.

🔗 Links:

# https://gdal.org/user/virtual_file_systems.html#vsizip-zip-archives
shape_url = "/vsizip/vsicurl/https://unosat-maps.web.cern.ch/MY/FL20191217MYS/FL20191217MYS_SHP.zip/ST1_20191215_WaterExtent_Johor_AOI2.shp"

This is a shapefile containing 🔷 polygons of the mapped water extent. Let’s put it into a DataPipe called zen3geo.datapipes.PyogrioReader (functional name: read_from_pyogrio).

dp_shapes = torchdata.datapipes.iter.IterableWrapper(iterable=[shape_url])
dp_pyogrio = dp_shapes.read_from_pyogrio()
dp_pyogrio

PyogrioReaderIterDataPipe

This will take care of loading the shapefile into a geopandas.GeoDataFrame object. Let’s take a look at the data table 📊 to see what attributes are inside.

it = iter(dp_pyogrio)
geodataframe = next(it)
geodataframe.dropna(axis="columns")
Water_Clas Sensor_ID Sensor_Dat Area_m2 Area_ha StaffID EventCode geometry
0 Flood Water Sentinel-1 2019-12-15 2.178154e+07 2178.154127 JT FL20191217MYS MULTIPOLYGON (((104.01442 1.72298, 104.01335 1...

Cool, and we can also visualize the polygons 🔷 on a 2D map. To align the coordinates with the 🛰️ Sentinel-1 image above, we’ll first use geopandas.GeoDataFrame.to_crs() to reproject the vector from 🌐 EPSG:9707 (WGS 84 + EGM96 height, latitude/longitude) to 🌏 EPSG:32648 (UTM Zone 48N).

print(f"Original bounds in EPSG:9707:\n{geodataframe.bounds}")
gdf = geodataframe.to_crs(crs="EPSG:32648")
print(f"New bounds in EPSG:32648:\n{gdf.bounds}")

Original bounds in EPSG:9707:
         minx      miny        maxx      maxy
0  103.844605  1.722825  104.187997  2.075693
New bounds in EPSG:32648:
            minx           miny           maxx           maxy
0  371483.597954  190459.328366  409683.869979  229473.310892

Plot it with geopandas.GeoDataFrame.plot(). This vector map 🗺️ should correspond to the zoomed in Sentinel-1 image plotted earlier above.

gdf.plot(figsize=(11.5, 9))

<Axes: >
_images/ffc378a80e3fd2c43f22a4d80b11481b2aff5f5191a0493426bd54c8e4ead3b6.png

Tip

Make sure to understand your raster and vector datasets well first! Open the files up in your favourite 🌐 Geographic Information System (GIS) tool, see how they actually look like spatially. Then you’ll have a better idea to decide on how to create your data pipeline. The zen3geo way puts you as the Master 🧙 in control.

1️⃣ Create a canvas to paint on 🎨#

In this section, we’ll work on converting the flood water 🌊 polygons above from a 🚩 vector to a 🌈 raster format, i.e. rasterization. This will be done in two steps 📶:

  1. Defining a blank canvas 🎞️

  2. Paint the polygons onto this blank canvas 🧑‍🎨

For this, we’ll be using tools from zen3geo.datapipes.datashader(). Let’s see how this can be done.

Blank canvas from template raster 🖼️#

A canvas represents a 2D area with a height and a width 📏. For us, we’ll be using a datashader.Canvas, which also defines the range of y-values (ymin to ymax) and x-values (xmin to xmax), essentially coordinates for every unit 🇾 height and 🇽 width.

Since we already have a Sentinel-1 🛰️ raster grid with defined height/width and y/x coordinates, let’s use it as a 📄 template to define our canvas. This is done via zen3geo.datapipes.XarrayCanvas (functional name: canvas_from_xarray).

dp_canvas = dp_decibel_canvas.canvas_from_xarray()
dp_canvas

XarrayCanvasIterDataPipe

Cool, and here’s a quick inspection 👀 of the canvas dimensions and metadata.

it = iter(dp_canvas)
canvas = next(it)
print(f"Canvas height: {canvas.plot_height}, width: {canvas.plot_width}")
print(f"Y-range: {canvas.y_range}")
print(f"X-range: {canvas.x_range}")
print(f"Coordinate reference system: {canvas.crs}")

Canvas height: 2743, width: 3538
Y-range: (95200.47858355992, 314640.4785835599)
X-range: (247773.19092160053, 530813.1909216006)
Coordinate reference system: EPSG:32648

This information should match the template Sentinel-1 dataarray 🏁.

print(f"Dimensions: {dict(dataarray.sizes)}")
print(f"Affine transform: {dataarray.rio.transform()}")
print(f"Bounding box: {dataarray.rio.bounds()}")
print(f"Coordinate reference system: {dataarray.rio.crs}")

Dimensions: {'band': 1, 'y': 2743, 'x': 3538}
Affine transform: | 80.00, 0.00, 247773.19|
| 0.00,-80.00, 314640.48|
| 0.00, 0.00, 1.00|
Bounding box: (247773.19092160053, 95200.47858355992, 530813.1909216006, 314640.4785835599)
Coordinate reference system: EPSG:32648

Rasterize vector polygons onto canvas 🖌️#

Now’s the time to paint or rasterize the vector geopandas.GeoDataFrame polygons 🔷 onto the blank datashader.Canvas! This would enable us to have a direct pixel-wise X -> Y mapping ↔️ between the Sentinel-1 image (X) and target flood label (Y).

The vector polygons can be rasterized or painted 🖌️ onto the template canvas using zen3geo.datapipes.DatashaderRasterizer (functional name: rasterize_with_datashader).

dp_datashader = dp_canvas.rasterize_with_datashader(vector_datapipe=dp_pyogrio)
dp_datashader

DatashaderRasterizerIterDataPipe

This will turn the vector geopandas.GeoDataFrame into a raster xarray.DataArray grid, with the spatial coordinates and bounds matching exactly with the template Sentinel-1 image 😎.

Note

Since we have just one Sentinel-1 🛰️ image and one raster 💧 flood mask, we have an easy 1:1 mapping. There are two other scenarios supported by zen3geo.datapipes.DatashaderRasterizer:

  1. N:1 - Many datashader.Canvas objects to one vector geopandas.GeoDataFrame. The single vector geodataframe will be broadcasted to match the length of the canvas list. This is useful for situations when you have a 🌐 ‘global’ vector database that you want to pair with multiple 🛰️ satellite images.

  2. N:N - Many datashader.Canvas objects to many vector geopandas.GeoDataFrame objects. In this case, the list of grids must ❗ have the same length as the list of vector geodataframes. E.g. if you have 5 grids, there must also be 5 vector files. This is so that a 1:1 pairing can be done, useful when each raster tile 🖽 has its own associated vector annotation.

See also

For more details on how rasterization of polygons work behind the scenes 🎦, check out Datashader’s documentation on:

2️⃣ Combine and conquer ⚔️#

So far, we’ve got two datapipes that should be 🧑‍🤝‍🧑 paired up in an X -> Y manner:

  1. The pre-processed Sentinel-1 🌈 raster image in dp_decibel_image

  2. The rasterized 💧 flood segmentation masks in dp_datashader

One way to get these two pieces in a Machine Learning ready chip format is via a stack, slice and split ™️ approach. Think of it like a sandwich 🥪, we first stack the bread 🍞 and lettuce 🥬, and then slice the pieces 🍕 through the layers once. Ok, that was a bad analogy, let’s just stick with tensors 🤪.

Stacking the raster layers 🥞#

Each of our 🌈 raster inputs are xarray.DataArray objects with the same spatial resolution and extent 🪟, so these can be stacked into an xarray.Dataset with multiple data variables. First, we’ll zip 🤐 the two datapipes together using torchdata.datapipes.iter.Zipper (functional name: zip)

dp_zip = dp_decibel_image.zip(dp_datashader)
dp_zip

ZipperIterDataPipe

This will result in a DataPipe where each item is a tuple of (X, Y) pairs 🧑‍🤝‍🧑. Just to illustrate what we’ve done so far, we can use torchdata.datapipes.utils.to_graph to visualize the data pipeline ⛓️.

torchdata.datapipes.utils.to_graph(dp=dp_zip)
_images/3d224fd099b1acd798574153e1a86197bbfcfb53acc34d8595d7fc0051bf4e27.svg

Next, let’s combine 🖇️ the two (X, Y) xarray.DataArray objects in the tuple into an xarray.Dataset using torchdata.datapipes.iter.Collator (functional name: collate). We’ll also ✂️ clip the dataset to a bounding box area where the target water mask has no 0 or NaN values.

def xr_collate_fn(image_and_mask: tuple) -> xr.Dataset:
    """
    Combine a pair of xarray.DataArray (image, mask) inputs into an
    xarray.Dataset with two data variables named 'image' and 'mask'.
    """
    # Turn 2 xr.DataArray objects into 1 xr.Dataset with multiple data vars
    image, mask = image_and_mask
    dataset: xr.Dataset = xr.merge(
        objects=[image.isel(band=0).rename("image"), mask.rename("mask")],
        join="override",
    )

    # Clip dataset to bounding box extent of where labels are
    mask_extent: tuple = mask.where(cond=mask == 1, drop=True).rio.bounds()
    clipped_dataset: xr.Dataset = dataset.rio.clip_box(*mask_extent)

    return clipped_dataset

dp_dataset = dp_zip.collate(collate_fn=xr_collate_fn)
dp_dataset

CollatorIterDataPipe

Double check to see that resulting xarray.Dataset’s image and mask looks ok 🙆‍♂️.

it = iter(dp_dataset)
dataset = next(it)

# Create subplot with VV image on the left and Water mask on the right
fig, axs = plt.subplots(ncols=2, figsize=(11.5, 4.5), sharey=True)
dataset.image.plot.imshow(ax=axs[0], cmap="Blues_r")
axs[0].set_title("Sentinel-1 VV channel")
dataset.mask.plot.imshow(ax=axs[1], cmap="Blues")
axs[1].set_title("Water mask")
plt.show()

/home/docs/checkouts/readthedocs.org/user_builds/zen3geo/envs/latest/lib/python3.11/site-packages/torch/utils/data/datapipes/iter/combining.py:297: UserWarning: Some child DataPipes are not exhausted when __iter__ is called. We are resetting the buffer and each child DataPipe will read from the start again.
  warnings.warn("Some child DataPipes are not exhausted when __iter__ is called. We are resetting "
_images/83edceb8ff0e10eeb39676748f583955d76fbfc234058d1fba2be4c5cdb71836.png

Slice into chips and turn into tensors 🗡️#

To cut 🔪 the xarray.Dataset into 128x128 sized chips, we’ll use zen3geo.datapipes.XbatcherSlicer (functional name: slice_with_xbatcher). Refer to Chipping and batching data if you need a 🧑‍🎓 refresher.

dp_xbatcher = dp_dataset.slice_with_xbatcher(input_dims={"y": 128, "x": 128})
dp_xbatcher

XbatcherSlicerIterDataPipe

Next step is to convert the 128x128 chips into a torch.Tensor via torchdata.datapipes.iter.Mapper (functional name: map). The 🛰️ Sentinel-1 image and 💧 water mask will be split out at this point too.

def dataset_to_tensors(chip: xr.Dataset) -> (torch.Tensor, torch.Tensor):
    """
    Converts an xarray.Dataset into to two torch.Tensor objects, the first one
    being the satellite image, and the second one being the target mask.
    """
    image: torch.Tensor = torch.as_tensor(chip.image.data)
    mask: torch.Tensor = torch.as_tensor(chip.mask.data.astype("uint8"))

    return image, mask

dp_map = dp_xbatcher.map(fn=dataset_to_tensors)
dp_map

MapperIterDataPipe

At this point, we could do some batching and collating, but we’ll point you again to Chipping and batching data to figure it out 😝. Let’s take a look at a graph of the complete data pipeline.

torchdata.datapipes.utils.to_graph(dp=dp_map)
_images/1376c62b77e9cc2f6b33ca84dd6185e837cf836f74d284bec443b86b3f468742.svg

Sweet, time for the final step ⏩.

Into a DataLoader 🏋️#

Pass the DataPipe into torch.utils.data.DataLoader 🤾!

dataloader = torch.utils.data.DataLoader(dataset=dp_map)
for i, batch in enumerate(dataloader):
    image, mask = batch
    print(f"Batch {i} - image: {image.shape}, mask: {mask.shape}")

/home/docs/checkouts/readthedocs.org/user_builds/zen3geo/envs/latest/lib/python3.11/site-packages/torch/utils/data/datapipes/iter/combining.py:297: UserWarning: Some child DataPipes are not exhausted when __iter__ is called. We are resetting the buffer and each child DataPipe will read from the start again.
  warnings.warn("Some child DataPipes are not exhausted when __iter__ is called. We are resetting "

Batch 0 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 1 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 2 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 3 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 4 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 5 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 6 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 7 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])
Batch 8 - image: torch.Size([1, 128, 128]), mask: torch.Size([1, 128, 128])

Now go train some flood water detection models 🌊🌊🌊

See also

To learn more about AI-based flood mapping with SAR, check out these resources: