Lab 04 assignment notebook 1 (20 pts)

Lab 04 assignment notebook 1 (20 pts)#

Notebook 1 of 2

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

Introduction#

This lab focuses on solidifying your understanding of fundamental raster concepts. You will…

inspect, read, and write raster data
understand raster attributes
perform quantitative analysis such as band math operations and thresholding
practice different raster visualization approaches

Instructions#

Please go through all cells in the notebook sequentially, making sure to complete all instructions
Some answers are left in as a guide–you still need to fill in the code to produce the output on your own!
Answers should be produced using code unless otherwise noted by “Written response”
Follow submission instructions closely, making sure you save and submit your notebook with all your outputs preserved

Part 0: Imports and filenames#

First, let’s import the libraries we’ll need.#

import geopandas as gpd
import rasterio as rio
import rioxarray as rxr
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar
from pathlib import Path
import os
import pyproj
import numpy as np

#May want to use interactive plotting for zoom/pan and live coordinate display
#%matplotlib widget
%matplotlib inline

Let’s make sure our files are where we think they are!#

If you haven’t downloaded the files already, please run the download demo notebook.

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

!ls -lh $imgdir

total 1.3G
-rw-r--r-- 1 eric eric 1.2M Jan 26 11:31 LC08_L2SP_046027_20180818_02_T1_rendered_preview.png
-rw-r--r-- 1 eric eric  81M Jan 26 11:27 LC08_L2SP_046027_20180818_20200831_02_T1_SR_B1.TIF
-rw-r--r-- 1 eric eric  81M Jan 26 11:28 LC08_L2SP_046027_20180818_20200831_02_T1_SR_B2.TIF
-rw-r--r-- 1 eric eric  82M Jan 26 11:28 LC08_L2SP_046027_20180818_20200831_02_T1_SR_B3.TIF
-rw-r--r-- 1 eric eric  83M Jan 26 11:29 LC08_L2SP_046027_20180818_20200831_02_T1_SR_B4.TIF
-rw-r--r-- 1 eric eric  92M Jan 26 11:29 LC08_L2SP_046027_20180818_20200831_02_T1_SR_B5.TIF
-rw-r--r-- 1 eric eric  90M Jan 26 11:30 LC08_L2SP_046027_20180818_20200831_02_T1_SR_B6.TIF
-rw-r--r-- 1 eric eric  84M Jan 26 11:31 LC08_L2SP_046027_20180818_20200831_02_T1_ST_B10.TIF
-rw-r--r-- 1 eric eric 115M Jan 27 19:39 LC08_L2SP_046027_20180818_20200831_02_T1_ST_B10_4326.TIF
-rw-r--r-- 1 eric eric 1.5M Jan 26 11:27 LC08_L2SP_046027_20181224_02_T1_rendered_preview.png
-rw-r--r-- 1 eric eric  91M Jan 26 11:23 LC08_L2SP_046027_20181224_20200829_02_T1_SR_B1.TIF
-rw-r--r-- 1 eric eric  91M Jan 26 11:24 LC08_L2SP_046027_20181224_20200829_02_T1_SR_B2.TIF
-rw-r--r-- 1 eric eric  92M Jan 26 11:24 LC08_L2SP_046027_20181224_20200829_02_T1_SR_B3.TIF
-rw-r--r-- 1 eric eric  91M Jan 26 11:25 LC08_L2SP_046027_20181224_20200829_02_T1_SR_B4.TIF
-rw-r--r-- 1 eric eric  96M Jan 26 11:26 LC08_L2SP_046027_20181224_20200829_02_T1_SR_B5.TIF
-rw-r--r-- 1 eric eric  90M Jan 26 11:26 LC08_L2SP_046027_20181224_20200829_02_T1_SR_B6.TIF
-rw-r--r-- 1 eric eric  77M Jan 26 11:27 LC08_L2SP_046027_20181224_20200829_02_T1_ST_B10.TIF

Let’s assign band filenames to variables#

#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

We’ve done the thermal infrared, shortwave infrared, and near infrared for you.

# B10 is the thermal infrared /surface temperture band
tir_august_fn = os.path.join(imgdir, august_id+'_ST_B10.TIF')
tir_december_fn = os.path.join(imgdir, december_id+'_ST_B10.TIF')

swir_august_fn = os.path.join(imgdir, august_id+'_SR_B6.TIF')
swir_december_fn = os.path.join(imgdir, december_id+'_SR_B6.TIF')

nir_august_fn = os.path.join(imgdir, august_id+'_SR_B5.TIF')
nir_december_fn = os.path.join(imgdir, december_id+'_SR_B5.TIF')

Please add variables for red, green, and blue each of the two images.#

If you’re unsure of which filename corresponds to each band, check the chart in the download demo notebook.

# STUDENT CODE HERE

Written response: Use `gdalinfo` to inspect the blue band from the december image, and report which CRS this image is in.#

# STUDENT CODE HERE

STUDENT WRITTEN RESPONSE HERE

Part 3: Raster transform (2 pts)#

Written response: Inspect the transform property. In your own words, what does this represent? How would it be used mathematically?#

transform = tir_august_da.rio.transform()
transform

Affine(59.99227997941328, 0.0, 473685.0,
       0.0, -59.99239736441966, 5373615.0)

STUDENT WRITTEN RESPONSE HERE

Use the transform to convert the raster’s corner pixel coordinates to projected coordinates.#

Hint: After you calculate the projected coordinates of the corners, try using .rio.bounds() to see if your answer make sense in the context of this raster’s bounds

upper_left_pixel = (0,0)
lower_left_pixel = (0,tir_august_da.rio.height)
upper_right_pixel = (tir_august_da.rio.width,0)
lower_right_pixel = (tir_august_da.rio.width,tir_august_da.rio.height)

print(f'Upper left pixel: {upper_left_pixel}')
print(f'Lower left pixel: {lower_left_pixel}')
print(f'Upper right pixel: {upper_right_pixel}')
print(f'Lower right pixel: {lower_right_pixel}')

Upper left pixel: (0, 0)
Lower left pixel: (0, 3946)
Upper right pixel: (3886, 0)
Lower right pixel: (3886, 3946)

# STUDENT CODE HERE

Upper left projected: (473685.0, 5373615.0)
Lower left projected: (473685.0, 5136885.0)
Upper right projected: (706815.0, 5373615.0)
Lower right projected: (706815.0, 5136885.0)

# STUDENT CODE HERE

(473685.0, 5136885.0, 706815.0, 5373615.0)

Use the corner projected coordinates to compute the total dimensions of the raster in km.#

Print your answer with an f string with a format similar to Total height: height km, total width: width km

# STUDENT CODE HERE

Part 5: Raster sampling (1 pts)#

Previously, rasterio users would have to use src.xy and src.index to convert between pixel and projected coordinates. They would either need to convert to pixel coordinates to index into the numpy array, or they would have to use src.sample.

Rioxarray’s (xarray’s) .sel() and .isel() functionalities can be really helpful for raster sampling as well. Since .isel() selects based integer array coordinates, you can enter in pixel coordinates. .sel() selects based on coordinate labels, so you can enter projected coordinates. When selecting projected coordinates which are often float values, we’ll use method='nearest' since we likely won’t know the exact coordinate discretization.

To select a point based on projected coordinates:

da.sel(x=x_coord,y=y_coord,method="nearest")

and for pixel coordinates:

da.isel(x=x_coord,y=y_coord)

Given the following points to sample, print the surface temperature value at that point. Note that the first two points are in the projected coordinates, the third point is in pixel coordinates, and the fourth point is in latitude and longitude.#

Hint: We’ve defined a Transformer for you to convert the 4th point’s latitude and longitude to a point in this CRS. See the sampling section of the rioxarray demo for more help.

glacier_proj_point = (595949,5190623)
sun_facing_proj_point = (600000,5195620)
shaded_pixel_point = (539,264)
parking_lot_latlon_point = (-121.7354,46.7856)

transform_4326_point_to_raster_crs = pyproj.Transformer.from_crs("EPSG:4326",tir_august_subset_da.rio.crs,always_xy=True)

# STUDENT CODE HERE

print(f'Glacier pixel temperature: {glacier_value:.2f}°C')
print(f'Sun-facing pixel temperature: {sun_facing_value:.2f}°C')
print(f'Shaded pixel temperature: {shaded_pixel_value:.2f}°C')
print(f'Paradise parking lot pixel temperature: {parking_lot_value:.2f}°C')

Glacier pixel temperature: 0.83°C
Sun-facing pixel temperature: 41.75°C
Shaded pixel temperature: 18.97°C
Paradise parking lot pixel temperature: 32.61°C

Challenge question: Plot the four points on top of the surface temperature raster and add a scalebar (GS: Attempt required / UG: +1 pts)#

Hint: You will likely have to transform the shaded pixel point and parking lot point before plotting

# STUDENT CODE HERE

../../_images/tir_band_points_of_interest.png

Written response: Do the temperatures seem reasonable for their respective locations?#

STUDENT WRITTEN RESPONSE HERE

Great work, you finished notebook 1! Move on to Notebook #2, but first…#

Save this notebook!
- You can also do an add, commit, and push at this stage
Shut down the kernel to free RAM
Proceed to Notebook #2!

Lab 04 assignment notebook 1 (20 pts)

Contents

Lab 04 assignment notebook 1 (20 pts)#

Introduction#

Instructions#

Part 0: Imports and filenames#

First, let’s import the libraries we’ll need.#

Let’s make sure our files are where we think they are!#

Let’s assign band filenames to variables#

Please add variables for red, green, and blue each of the two images.#

Written response: Use gdalinfo to inspect the blue band from the december image, and report which CRS this image is in.#

Part 1: Raster basics (2 pts)#

Before we switch to rioxarray for the remainder of the lab, let’s briefly open and read in the near infrared band from the August 2018 image with rasterio.#

Written response: Based on this inspection of the metadata, how do you think missing values are represented?#

Unmasked…#

Masked…#

Written response: Compare the amount of memory the unmasked array takes up in memory vs how much memory the respective file takes up on the disk (see !ls -lh $imgdir output above). Why might these be different?#

Written response: What differences did you notice between the unmasked and masked ouputs and plots? Do the dark purple regions in the unmasked array reflect real values?#

Now we’ll do the same thing, but using rioxarray…#

Written response: Based on your observations, describe how rasterio and rioxarray handle masking differently. Based on only this insight, give an example of a case in which you might prefer to use rasterio over rioxarray.#

Inspect and print the following raster metadata using the xarray rio accessor… shape, height, width, nodata value, encoded nodata value, number of bands, CRS, resolution, extent, and transform. Try to print things neatly if you can using f strings.#

What is the value of the right bound?#

Now we’ll open the red band for August with rioxarray, using masked=True and overview_level=2.#

Written response: What did setting overview_level=2 do?#

Make a figure with two subplots. On the left, plot the raster using the underlying numpy array. On the right, plot the raster using xarray DataArray’s built in .plot.imshow() functionality. Use a gray colormap with the same vmin and vmax for both, and please include a title on each.#

Written response: Why does the Landsat-8 image appear rotated relative to the axes?#

Written response: What is different about the two axes? We know the dimensions of the two arrays are the same, so what does this tell us about xarray DataArray coordinates? Feel free to inspect the xarray DataArray for clues.#

Part 2: Converting to geophysical variables (3 pts)#

First, let’s open the thermal infrared band with masked=False and overview=0…#

Make a figure with a plot of the thermal infrared image on the left, and its corresponding histogram on the right.#

Written response: Why do you see a spike at 0 in the histogram?#

Try making the previous figure again, this time loading with masked=True#

These values still don’t look very interpretable. Let’s convert them to degrees C.#

Calculate the min and max surface temperature values, and the 2nd and 98th percentile of surface temperature values.#

With your updated surface temperature values, please make a final figure with pixel values on the left and the histogram on the right. Try using the inferno colormap, and set vmin=0 and vmax=40. Please also draw vertical lines for the 2nd and 98th percentile on the histogram.#

Written response: Do the values make sense for August? Based on your image, where are the temperatures generally warmest? Coldest?#

Part 3: Raster transform (2 pts)#

Written response: Inspect the transform property. In your own words, what does this represent? How would it be used mathematically?#

Use the transform to convert the raster’s corner pixel coordinates to projected coordinates.#

Use the corner projected coordinates to compute the total dimensions of the raster in km.#

Part 4: Raster windows (1 pts)#

Let’s try extracting a window around Mt. Rainier using both methods#

Use .sel() and .isel() to select a window around Mt. Rainier. Convert to surface temperature in degrees C, and then plot these next to each other.#

Part 5: Raster sampling (1 pts)#

Given the following points to sample, print the surface temperature value at that point. Note that the first two points are in the projected coordinates, the third point is in pixel coordinates, and the fourth point is in latitude and longitude.#

Challenge question: Plot the four points on top of the surface temperature raster and add a scalebar (GS: Attempt required / UG: +1 pts)#

Written response: Do the temperatures seem reasonable for their respective locations?#

Great work, you finished notebook 1! Move on to Notebook #2, but first…#

Written response: Use `gdalinfo` to inspect the blue band from the december image, and report which CRS this image is in.#

Written response: Compare the amount of memory the unmasked array takes up in memory vs how much memory the respective file takes up on the disk (see `!ls -lh $imgdir` output above). Why might these be different?#

Inspect and print the following raster metadata using the xarray `rio` accessor… shape, height, width, nodata value, encoded nodata value, number of bands, CRS, resolution, extent, and transform. Try to print things neatly if you can using f strings.#

Now we’ll open the red band for August with rioxarray, using `masked=True` and `overview_level=2`.#

Written response: What did setting `overview_level=2` do?#

Make a figure with two subplots. On the left, plot the raster using the underlying numpy array. On the right, plot the raster using xarray DataArray’s built in `.plot.imshow()` functionality. Use a gray colormap with the same vmin and vmax for both, and please include a title on each.#

First, let’s open the thermal infrared band with `masked=False` and `overview=0`…#

Try making the previous figure again, this time loading with `masked=True`#

With your updated surface temperature values, please make a final figure with pixel values on the left and the histogram on the right. Try using the `inferno` colormap, and set `vmin=0` and `vmax=40`. Please also draw vertical lines for the 2nd and 98th percentile on the histogram.#

Use `.sel()` and `.isel()` to select a window around Mt. Rainier. Convert to surface temperature in degrees C, and then plot these next to each other.#