Raster packages demo: rioxarray (3/3)

Raster packages demo: rioxarray (3/3)#

UW Geospatial Data Analysis
CEE467/CEWA567
David Shean, Eric Gagliano, Quinn Brencher

Introduction#

Working with raster data is a essential part of geospatial analysis. There are several python packages available for raster processing, and we’ll introduce three key tools that build on top of each other…

GDAL (Geospatial Data Abstraction Library): The low-level foundation that powers most geospatial software. While powerful, it can be complex to use directly.
rasterio: A more Pythonic interface to GDAL that provides efficient access to raster data using NumPy arrays.
rioxarray: Higher-level package that combines the power of rasterio with xarray’s labeled dimensions and advanced capabilities for handling multi-dimensional data.

Each package plays an important role in the python geospatial ecosystem, so we’ll briefly introduce the tools one at a time to practice some fundamentals and gain some raster intuition.

The lower level GDAL and rasterio are very well-supported, and there are indeed use cases for when you might prefer interacting with these lower level tools. Ultimately, we’ll focus on rioxarray for the rest of the quarter due to its intuitive handling of multi-dimensional data (e.g. raster time series) and dask integration for scalability.

import os
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar
import rioxarray as rxr
import xarray as xr
import pyproj

imgdir = f'{Path.home()}/gda_demo_data/LS8_data'

#Pre-identified cloud-free Image IDs used for the lab
august_id = 'LC08_L2SP_046027_20180818_20200831_02_T1' # August 2018
december_id = 'LC08_L2SP_046027_20181224_20200829_02_T1' # # December 2018

# B2 is the blue band
blue_august_fn = os.path.join(imgdir, august_id+'_SR_B2.TIF')
blue_december_fn = os.path.join(imgdir, december_id+'_SR_B2.TIF')

# define a filename for the reprojected image
blue_august_4326_fn = os.path.join(imgdir, august_id+'_SR_B2_4326.TIF')

print(blue_august_fn)

/home/jovyan/gda_demo_data/LS8_data/LC08_L2SP_046027_20180818_20200831_02_T1_SR_B2.TIF

rioxarray#

Use the rioxarray open_rasterio function.

This is an xarray DataArray. The DataArray object contains our raster as a numpy array, but also contains the geographic coordinates and metadata.

blue_august_da.rio.crs

CRS.from_wkt('PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]')

blue_august_da.rio.nodata

np.float32(nan)

blue_august_da.rio.encoded_nodata # what's the difference?

np.float32(0.0)

#These are the standard scale and offset values
#SR 0.0000275 + -0.2
sr_scale = 0.0000275
sr_offset = -0.2
#ST 0.00341802 + 149.0
st_scale = 0.00341802
st_offset = 149.0

blue_august_da.values

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

blue_august_da.plot.imshow(cmap='Blues', vmin=0, vmax=0.4)

<matplotlib.image.AxesImage at 0x7f0580a3e6d0>

../../_images/ac46fd499ccbb7d836f70a1697dd9382e157f2a0d71df14cc8f1907eea74856c.png

# try adding a scalebar
f,ax=plt.subplots()
blue_august_da.plot.imshow(ax=ax,cmap='Blues', vmin=0, vmax=0.4)
ax.set_aspect('equal')
ax.add_artist(ScaleBar(1))

<matplotlib_scalebar.scalebar.ScaleBar at 0x7f058099ad50>

../../_images/9f99968c801ce919a0d67f129a6ded83bbbcf4143fda17aab9de99754606f552.png

f,ax=plt.subplots()
ax.hist(blue_august_da.values.ravel(),bins=100);

../../_images/8c299a87b1e42ec244951e8d95acc2bbc59b085b9e7729547e3087be01a9c80f.png

Sample some points in the raster#

lumen_field_proj_point = (550256, 5271504)
husky_stadium_pixel_point = (2623,3198)
lake_washington_latlon_point = (-122.25, 47.6)

transform_4326_point_to_raster_crs = pyproj.Transformer.from_crs("EPSG:4326",blue_august_da.rio.crs,always_xy=True)
parking_lot_latlon_point = transform_4326_point_to_raster_crs.transform(lake_washington_latlon_point[0],lake_washington_latlon_point[1])

lumen_field_value = float(blue_august_da.sel(x=lumen_field_proj_point[0],y=lumen_field_proj_point[1],method='nearest'))
husky_stadium_value = float(blue_august_da.isel(x=husky_stadium_pixel_point[0],y=husky_stadium_pixel_point[1]))
lake_washington_value = float(blue_august_da.sel(x=parking_lot_latlon_point[0],y=parking_lot_latlon_point[1],method='nearest'))

print(f"Lumen Field: {lumen_field_value}")
print(f"Husky Stadium: {husky_stadium_value}")
print(f"Lake Washington: {lake_washington_value}")

Lumen Field: 0.08088500797748566
Husky Stadium: 0.11058498919010162
Lake Washington: 0.0016299933195114136

Try reprojecting#

Pretty easy!

# reproject
blue_august_da.rio.reproject('EPSG:4326').rio.to_raster(blue_august_4326_fn)