Rasterizing vectors & vectorizing rasters

• Sign up to the DEA Sandbox to run this notebook interactively from a browser

• Compatibility: Notebook currently compatible with both the NCI and DEA Sandbox environments

• Products used: ga_ls_wo_fq_cyear_3

Background

Many remote sensing and/or geospatial workflows require converting between vector data (e.g. shapefiles) and raster data (e.g. pixel-based data like that in an xarray.DataArray). For example, we may need to use a shapefile as a mask to limit the analysis extent of a raster, or have raster data that we want to convert into vector data to allow for easy geometry operations.

Description

This notebook demonstrates the use of the DEA function xr_rasterize and xr_vectorize from dea_tools.spatial.

The first section loads in DEA Water Observations Statistics data from Digital Earth Australia, and vectorises the pixel-based xarray.DataArray object into a vector-based geopandas.GeoDataFrame object containing persistent water-bodies as polygons. We then export the GeoDataframe as a vector file.

The second section rasterises the vector data we created in the first section back into an xarray.DataArray, and exports the results as a GeoTIFF.

Getting started

To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.

Load packages

import datacube
import matplotlib.pyplot as plt

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.spatial import xr_vectorize, xr_rasterize

Connect to the datacube

dc = datacube.Datacube(app='Rasterize_vectorize')

Load DEA Water Observations data from the datacube

We will load in an annual summary from the DEA Water Observations Statistics product to provide us with some data to work with. The query below will load the 2000 annual summary of DEA Water Observations product for a region around the Menindee Lakes.

# Create a query object
query = {
    'x': (142.1, 142.80),
    'y': (-32.1, -32.6),
    'time': ('2000')
}

# Load WoFS through the datacube
ds = dc.load(product='ga_ls_wo_fq_cyear_3',
             **query)
ds

<xarray.Dataset>
Dimensions:      (time: 1, y: 2037, x: 2324)
Coordinates:
  * time         (time) datetime64[ns] 2000-07-01T23:59:59.999999
  * y            (y) float64 -3.536e+06 -3.536e+06 ... -3.597e+06 -3.597e+06
  * x            (x) float64 9.396e+05 9.396e+05 ... 1.009e+06 1.009e+06
    spatial_ref  int32 3577
Data variables:
    count_wet    (time, y, x) int16 0 0 0 0 0 0 0 0 0 1 ... 0 0 0 0 0 0 0 0 0 0
    count_clear  (time, y, x) int16 29 29 29 29 29 29 29 ... 14 14 14 14 14 14
    frequency    (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Attributes:
    crs:           EPSG:3577
    grid_mapping:  spatial_ref

xarray.Dataset

• Dimensions:
  • time: 1
  • y: 2037
  • x: 2324
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2000-07-01T23:59:59.999999

    units :

    seconds since 1970-01-01 00:00:00

    array(['2000-07-01T23:59:59.999999000'], dtype='datetime64[ns]')

  • y
    (y)
    float64
    -3.536e+06 ... -3.597e+06

    units :

    metre

    resolution :

    -30.0

    crs :

    EPSG:3577

    array([-3535845., -3535875., -3535905., ..., -3596865., -3596895., -3596925.])

  • x
    (x)
    float64
    9.396e+05 9.396e+05 ... 1.009e+06

    units :

    metre

    resolution :

    30.0

    crs :

    EPSG:3577

    array([ 939615.,  939645.,  939675., ..., 1009245., 1009275., 1009305.])

  • spatial_ref
    ()
    int32
    3577

    spatial_ref :

    PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]

    grid_mapping_name :

    albers_conical_equal_area

    array(3577, dtype=int32)

• Data variables: (3)
  • count_wet
    (time, y, x)
    int16
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    1

    nodata :

    -999

    crs :

    EPSG:3577

    grid_mapping :

    spatial_ref

    array([[[0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            ...,
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0]]], dtype=int16)

  • count_clear
    (time, y, x)
    int16
    29 29 29 29 29 ... 14 14 14 14 14

    units :

    1

    nodata :

    -999

    crs :

    EPSG:3577

    grid_mapping :

    spatial_ref

    array([[[29, 29, 29, ..., 14, 14, 14],
            [29, 29, 29, ..., 14, 14, 14],
            [29, 29, 29, ..., 14, 14, 14],
            ...,
            [24, 24, 24, ..., 14, 14, 14],
            [24, 24, 24, ..., 14, 14, 14],
            [24, 24, 24, ..., 14, 14, 14]]], dtype=int16)

  • frequency
    (time, y, x)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    1

    nodata :

    nan

    crs :

    EPSG:3577

    grid_mapping :

    spatial_ref

    array([[[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]]], dtype=float32)

• Attributes: (2)

  crs :

  EPSG:3577

  grid_mapping :

  spatial_ref

Plot the DEA Water Observations data

Let’s plot the DEA Water Observations data to get an idea of the objects we will be transforming. In the code below, we first select the pixels where the satellite has observed water at least 25% of the year, this is so we can isolate the more persistent water bodies and reduce some of the noise before we vectorize the raster.

# Select pixels that are classified as water > 25 % of the time
water_bodies = ds.frequency > 0.25

# Plot the data
water_bodies.plot(size=5)

<matplotlib.collections.QuadMesh at 0x7fa3b9a7ba30>

Vectorizing an xarray.DataArray

To convert our xarray.DataArray object into a vector based geopandas.GeoDataFrame, we can use the DEA function xr_vectorize from dea_tools.spatial. This tool is based on the rasterio.features.shape function, and can accept any of the arguments in rasterio.features.shape using the same syntax.

In the cell below, we use the argument mask=water_bodies.values==1 to indicate we only want to convert the values in the xarray object that are equal to 1.

gdf = xr_vectorize(da=water_bodies,
                   mask=water_bodies.values == 1)

Plot our vectorised raster

gdf.plot(figsize=(6, 6))

<AxesSubplot:>

Export as vector file

The xr_vectorize function also allows us to very easily export the GeoDataFrame as a vector file for use in other applications using the output_path parameter. For this example we will export our data as a GeoJSON vector file; use the extension .shp to export an ESRI Shapefile instead.

gdf = xr_vectorize(da=water_bodies,
                   mask=water_bodies.values == 1,
                   output_path='test.geojson')

Exporting vector data to test.geojson

---

Rasterizing a vector file

Using the xr_rasterize function from `dea_tools.spatial <../Tools/dea_tools/spatial.py>`__, we can “rasterize” our geopandas.GeoDataFrame back into a xarray.DataArray. This function uses an existing xarray.DataArray object as a template for converting our vector into a raster object (the template provides the size, CRS, dimensions, and attributes of the output array).

As we already have the GeoDataFrame loaded we don’t need to read in our vector file, but if we wanted to read in a vector file first we could use `gpd.read_file() <https://geopandas.org/en/stable/docs/reference/api/geopandas.read_file.html>`__.

Note: The xr_rasterize function is based on the rasterio function `rasterio.features.rasterize <https://rasterio.readthedocs.io/en/stable/api/rasterio.features.html>`__, and can accept any arguments supported by rasterio.features.rasterize.

water_bodies_again = xr_rasterize(gdf=gdf, da=water_bodies)

We can plot our rasterised water_bodies_again data to verify it looks identical to our original water_bodies data we loaded from the datacube:

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
water_bodies.plot(ax=axes[0])
water_bodies_again.plot(ax=axes[1])
axes[0].set_title('Original waterbodies')
axes[1].set_title('Rasterized waterbodies');

Rasterizing attribute values

In the previous examples, we rasterized our vector data into a simple boolean raster containing 1s and 0s. We can however also use a specific attribute column in our vector data as values that will be assigned to each of our raster pixels.

For example, we can modify our vectorised data to give each polygon a unique number in the ‘attribute’ column:

gdf['attribute'] = gdf.index
gdf

	attribute	geometry
0	0	POLYGON ((942810.000 -3535830.000, 942810.000 ...
1	1	POLYGON ((942930.000 -3535830.000, 942930.000 ...
2	2	POLYGON ((943140.000 -3535860.000, 943200.000 ...
3	3	POLYGON ((1000410.000 -3535830.000, 1000440.00...
4	4	POLYGON ((950550.000 -3536010.000, 950550.000 ...
...	...	...
1293	1293	POLYGON ((950550.000 -3596790.000, 950550.000 ...
1294	1294	POLYGON ((946980.000 -3592680.000, 947190.000 ...
1295	1295	POLYGON ((942810.000 -3596670.000, 942840.000 ...
1296	1296	POLYGON ((942780.000 -3596910.000, 942780.000 ...
1297	1297	POLYGON ((964620.000 -3596730.000, 964650.000 ...

1298 rows × 2 columns

Now we can rasterize our polygons using the unique values from the ‘attribute’ column:

water_bodies_attribute = xr_rasterize(gdf=gdf,
                                      da=water_bodies,
                                      attribute_col='attribute')

When we plot this raster data, we can see our rasterized waterbodies now contain pixel values that match the unique values in the ‘attribute’ column of our vector data:

water_bodies_attribute.plot(cmap='nipy_spectral')

<matplotlib.collections.QuadMesh at 0x7fa3addb8520>

Export rasterized data as a GeoTIFF

xr_rasterize also allows us to export our results as a GeoTIFF raster file using the output_path parameter:

water_bodies_again = xr_rasterize(gdf=gdf,
                                  da=water_bodies,
                                  output_path='test.tif')

Exporting raster data to test.tif

---

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Discord chat or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on GitHub.

Last modified: December 2023

Compatible datacube version:

print(datacube.__version__)

1.8.6

Tags

Tags: sandbox compatible, NCI compatible, DEA Water Observations, rasterize, rasterio.features.shapes, rasterio.features.rasterize, GeoPandas, vectorize, GeoTIFF, vector, xr_rasterize, xr_vectorize