Skip to article frontmatterSkip to article content

Download the sample datasets for this tutorial

For ease of access during the hackweek, sample files are available for download for running the command in the cell below. (See the bonus notebook “thermal-ir-data-download.ipynb” for more details about data access methods)

# S3 base URL for tutorial data
S3_BASE_URL = "https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/"

Import the packages we’ll need for this tutorial

# Import some general-purpose packages for handling different data structures
import numpy as np # for working with n-D arrays
import pandas as pd # for reading our csv data file and working with tabular data

# Import matplotlib which we'll use for plotting images and graphs
import matplotlib.pyplot as plt

# Import these packages for working with raster data
import xarray as xr # xarray lets us work with n-D arrays and labeled data, such as NetCDF files
import rioxarray # rioxarray provides capabilities of the rasterio package to xarray, letting us easily work with files such as GeoTIFFs
import fsspec

# Import some packages for working with the SnowEx SQL database
from snowexsql.db import get_db # Import the connection function from the snowexsql library
from  snowexsql.data import SiteData # Import the table classes from our data module which is where our ORM classes are defined 
from datetime import datetime # Import some tools to build dates 
from snowexsql.conversions import query_to_geopandas # Import a useful function for plotting and saving queries! See https://snowexsql.readthedocs.io/en/latest/snowexsql.html#module-snowexsql.conversions

Part 1: Comparing airborne IR imagery with ground-truth observations

Airborne IR imagery

The Naval Postgraduate School Twin Otter aircraft carried the UW APL thermal infrared imager and SWESARR instrument over Grand Mesa for SnowEx 2020. (Photo by Chris Chickadel)

Load IR image mosaic geotiff file

airborne_ir = rioxarray.open_rasterio(S3_BASE_URL + 'SNOWEX2020_IR_PLANE_2020Feb08_mosaicked_2020-02-08T181915.tif')

Inspect the contents of the file we just opened

airborne_ir
Loading...
airborne_ir.rio.crs # original CRS
CRS.from_wkt('PROJCS["WGS 84 / UTM zone 12N",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",-111],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","32612"]]')

Under the dataarray’s attributes we can see that we have a coordinate reference system already defined (crs) as EPSG:32612. We can also find this through da.rio.crs.

However, we would like to reproject this into the common projection used by datasets on the SnowEx SQL database: EPSG:26912. We can do that using rioxarray’s reproject method (see an example here).

NOTE: Reprojections typically require a lot of processing time, so expect this cell to run for up to 1 minute.

airborne_ir = airborne_ir.rio.reproject('EPSG:26912') # overwrite itself with new reprojected data array
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/SNOWEX2020_IR_PLANE_2020Feb08_mosaicked_2020-02-08T181915.tif.msk: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/SNOWEX2020_IR_PLANE_2020Feb08_mosaicked_2020-02-08T181915.tif.MSK: 403

airborne_ir.rio.crs # new CRS
CRS.from_wkt('PROJCS["NAD83 / UTM zone 12N",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-111],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","26912"]]')

Next, the filename shows us when this imagery was taken in UTC time, “2020-02-08T181915”

We can create a pandas timestamp variable in local time for comparison with other datasets:

# Create a pandas timestamp, subtract 7 hours from UTC time to get local time (MST, UTC-7)
airborne_ir_timestamp = pd.Timestamp(2020,2,8,18,19,15) - pd.Timedelta(hours=7)
What color scale should we use for temperature?

Common advice you may have heard is to avoid using rainbow color scales. Luckily matplotlib gives us lots of options to choose from.

When representing images of temperature, sometimes we want to pick colors that intuitively suggest temperature, such as the “magma” colorbar below. Other times we might be interested in both magnitude and sign, such as temperatures above or below melting point, in which case we could use something like “RdBu_r” below. I often pick simple greyscale, though less visually interesting, it is sometimes easier to pick out details in a continuous color scale like they “Greys” scale below. Make sure to include a labeled colorbar so your plot can be correctly interpreted!

Some matplotlib color scale options (top to bottom): “magma”, “RdBu_r”, “Greys”

Plot the airborne TIR image.

  • What does our chosen colorscale tell us about temperatures here?

  • What can we see?

fig, ax = plt.subplots(figsize=(20,10)) # create a new matplotlib figure, set the figure size
ax.set_aspect('equal') # set the aspect ratio to "equal"

# plot the airborne infrared image
airborne_ir.plot(cmap='magma', vmin=-10, vmax=10,ax=ax, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})

# set axes labels
ax.set_xlabel('Eastings UTM 12N (m)')
ax.set_ylabel('Northings UTM 12N (m)')

ax.set_title('Airborne TIR imagery of Grand Mesa, CO (Feb. 8, 2020)');
<Figure size 2000x1000 with 2 Axes>

Bonus activity:

To help with our image interpretation, we can load visible imagery taken concurrently from the UW-APL airborne instrument. (Note: this example image is a single band black and white image, though we also have full RGB images available through NSIDC)

airborne_vis = rioxarray.open_rasterio(S3_BASE_URL + 'SNOWEX2020_EO_PLANE_2020Feb08_mosaicked_2020-02-08T181915.tif')

# note that the filename is identical with the same timestamp, but is labeled "EO" (electro-optical) rather than "IR" (infrared)
# Also reproject the airborne visible imagery into EPSG:26912
airborne_vis = airborne_vis.rio.reproject('EPSG:26912') # overwrite itself with new reprojected data array
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/SNOWEX2020_EO_PLANE_2020Feb08_mosaicked_2020-02-08T181915.tif.msk: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/SNOWEX2020_EO_PLANE_2020Feb08_mosaicked_2020-02-08T181915.tif.MSK: 403

airborne_vis
Loading...

Plot the visible and infrared images side by side. This time, we will change the x and y axes limits (set_xlim, and set_ylim) to zoom in closer.

fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(20,10), tight_layout=True)

# Plot the IR imagery
airborne_ir.plot(ax=axs[0], 
                 cmap='magma', 
                 vmin=-10, vmax=10, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[0].set_title('Airborne IR')

# Plot the visible imagery
airborne_vis.plot(ax=axs[1], 
                  cmap='Greys_r',
                  cbar_kwargs={'label': 'DN'})
axs[1].set_title('Airborne Vis')

# for each subplot, do the following:
for ax in axs: 
    ax.set_aspect('equal') # set the aspect ratio to "equal"
    
    # give each axis a label
    ax.set_xlabel('Eastings UTM 12N (m)')
    ax.set_ylabel('Northings UTM 12N (m)')
    
    # set the axes limits, units in meters UTM Zone 12N (I chose these values by just looking at the plot above)
    ax.set_xlim((735000, 760000)) # x axis limits
    ax.set_ylim((4320000, 4325000)) # y axis limits
<Figure size 2000x1000 with 4 Axes>

Ground-based temperature observations

To provide a source of “ground truth” for the airborne and satellite thermal infrared images during the SnowEx 2020 Grand Mesa campaign, we can use ground-based snow surface temperature measurements. On February 5, 2020, we installed a thermal infrared radiometer pointing at the snow surface at snow pit #2S10 (left), and buried temperature sensors beneath the snow surface (right). These logged observations at 5-minute intervals until we removed the instrumentation a week later on February 12.

Snow temperature sensor setup at snow pit 2S10: (left) tripod-mounted thermal ifrared radiometer to measure snow surface, (right) temperature probes to be buried beneath the snow surface.
(Photos by Steven Pestana)

What are some of the differences we might expect to see between the ground-based surface temperature data and the thermal IR images?

  • Emissivity differences?

    • Our ground-based radiometer was looking at 45 deg off nadir, versus nadir ASTER versus variable view angle airborne (snow emissivity changes off-nadir)

  • Different TIR bandwidths (broad versus narrow)?

    • Ground-based radiometer: 8-14 μm, Airborne IR cameras: 8-14 μm, ASTER band 14: 10.95-11.65 µm

  • Different atmospheric path lengths?

    • From <1 meter, to 1 km, to entire atmospheric column (~100 km)

  • “Point” versus area, and geolocation accuracy

    • The ground-based radiometer is measuring temperature for a spot (not really a single point) maybe ~1m in diameter. The airborne and ASTER imagers have spatial resolutions of 5m and 90m respectively. How confident are we in the geolocation of individual pixels in the imagery?

Where is snow pit 2S10?

We can find this information through a query to the SnowEx SQL database. First, set up the connection:

# This is what you will use for all of hackweek to access the db
db_name = 'snow:hackweek@db.snowexdata.org/snowex'
# Using the function get_db, we receive 2 ways to interact with the database
engine, session = get_db(db_name)

Then, query SiteData using filter_by to find the entry with the site ID that we want (2S10). Preview the resulting geodataframe.

# Form the query to receive site_id='2S10' from the sites table
qry = session.query(SiteData).filter_by(site_id='2S10')

# Convert the record received into a geopandas dataframe
siteData_df = query_to_geopandas(qry, engine)

# Preview the resulting geopandas dataframe
siteData_df
Loading...

Inspect of the geodatframe’s metadata

# What is the coordinate reference system used here?
siteData_df.crs
<Projected CRS: EPSG:26912> Name: NAD83 / UTM zone 12N Axis Info [cartesian]: - E[east]: Easting (metre) - N[north]: Northing (metre) Area of Use: - name: North America - between 114°W and 108°W - onshore and offshore. Canada - Alberta; Northwest Territories; Nunavut; Saskatchewan. United States (USA) - Arizona; Colorado; Idaho; Montana; New Mexico; Utah; Wyoming. - bounds: (-114.0, 31.33, -108.0, 84.0) Coordinate Operation: - name: UTM zone 12N - method: Transverse Mercator Datum: North American Datum 1983 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
# Preview the geometry of this geodataframe, we should see that it is a POINT
siteData_df.geometry
0 POINT (743076 4322689) 1 POINT (743076 4322689) 2 POINT (743076 4322689) 3 POINT (743076 4322689) 4 POINT (743076 4322689) Name: geom, dtype: geometry

We can now plot our snow pit site from this geodataframe on top of the airborne IR image. (See more tips about plotting geodataframes here)

fig, ax = plt.subplots(figsize=(20,5)) # create a new matplotlib figure, set the figure size
ax.set_aspect('equal') # set the aspect ratio to "equal"

# plot the airborne infrared image
airborne_ir.plot(cmap='magma', vmin=-10, vmax=10, ax=ax, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})

# plot the location of the snow pit of interest
siteData_df.plot(ax=ax, color='r', marker='x')

# set axes labels
ax.set_xlabel('Eastings UTM 12N (m)')
ax.set_ylabel('Northings UTM 12N (m)')

# set the axes limits, units in meters UTM Zone 12N (I chose these values by just looking at the plot above)
ax.set_xlim((735000, 760000)) # x axis limits
ax.set_ylim((4320000, 4325000)) # y axis limits

# set plot title
ax.set_title('Location of Snow Pit 2S10\nwith Airborne TIR imagery of Grand Mesa, CO (Feb. 8, 2020)');
<Figure size 2000x500 with 2 Axes>

Change the x and y axes limits (set_xlim, and set_ylim) to zoom in to our point of interest. In this case we can use df.geometry.total_bounds to get the x and y values that define the area our geometry takes up. (In this case we have a point so it will return just the point’s location, but this would work if we had a polygon as well)

fig, ax = plt.subplots(figsize=(10,10)) # create a new matplotlib figure, set the figure size
ax.set_aspect('equal') # set the aspect ratio to "equal"

# plot the airborne infrared image
airborne_ir.plot(cmap='magma', vmin=-10, vmax=10, ax=ax, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})

# plot the location of the snow pit of interest
siteData_df.plot(ax=ax, color='r', marker='x')

# set axes limits
xmin, ymin, xmax, ymax = siteData_df.geometry.total_bounds # get the "total bounds" for our geometry
ax.set_xlim((xmin-1000, xmax+1000)) # x axis limits to +/- 1 km from our point's "total bounds"
ax.set_ylim((ymin-1000, ymax+1000)) # y axis limits to +/- 1 km from our point's "total bounds"

# set axes labels
ax.set_xlabel('Eastings UTM 12N (m)')
ax.set_ylabel('Northings UTM 12N (m)')

# set plot title
ax.set_title('Location of Snow Pit 2S10\nAirborne TIR imagery of Grand Mesa, CO (Feb. 8, 2020)');
<Figure size 1000x1000 with 2 Axes>

Import the snow temperature timeseries dataset

This data is available through NSIDC, but we have already downloaded a local copy for this tutorial. (See the bonus notebook thermal-ir-data-download.ipynb for more details about data access methods)

The raw data file doesn’t include the column names, so we need to set the column headers following the dataset’s README file.

!cat data/snow-temperature-README.txt
Description:
	This folder contains data files from two Campbell Scientific CR10X dataloggers which were installed at two sites on Grand Mesa to measure snow temperature profiles as part of the SnowEx 2020 field campaign.
	
	"Grand Mesa 1" (GM1) measured snow surface temperature with an Apogee SI-111 thermal infrared radiometer mounted on a tripod at a height of about 127 cm above the snow surface. The radiometer pointed 45 degrees off-nadir at the snow surface just to the west of the tripod. Five Campbell Scientific temperature probes were mounted vertically along a wooden stake with zip ties, at depths of 5, 10, 15, 20, and 30 cm below the snow surface. (CR10X_GM1_final_storage_1.dat) 
	
	"Grand Mesa 2" (GM2) measured snow temperatures with five temperature probes mounted vertically on a wooden stake with zip ties at depths of 5, 10, 15, 20, and 30 cm below the snow surface. This location was just outside the fence perimeter of the Mesa West met station (northwest of the met station tower). Snow surface temperature is measured at the Mesa West station by its own set of thermal infrared radiometers, which can provide snow surface temperature adjacent to this vertical profile of snow temperature. The Mesa West station should also be able to provide similar snow temperatures below the surface with its own string of temperature sensors.

Date/times: 
	GM1: 2-5-2020 11:12 AM MST (UTC−07:00) - 2-12-2020 3:36 PM MST (UTC−07:00)
	GM2: 2-5-2020 12:35 PM MST (UTC−07:00) - 2-12-2020 2:36 PM MST (UTC−07:00)

Locations: 
	GM1 was located at Pit ID 2S10
	GM2 was located at the Mesa West met station 


File formats:
*.dat - comma-separated values from the Campbell Scientific CR10X Dataloggers
*.xlsx - MS Excel worksheet

Columns:
	GM1:
		'table', 'year', 'doy', 'time',
		'rad_avg', 'rad_max', 'rad_min', 'rad_std',
		'sb_avg', 'sb_max', 'sb_min', 'sb_std',
		'temp1_avg', 'temp1_max', 'temp1_min', 'temp1_std',
		'temp2_avg', 'temp2_max', 'temp2_min', 'temp2_std',
		'temp3_avg', 'temp3_max', 'temp3_min', 'temp3_std',
		'temp4_avg', 'temp4_max', 'temp4_min', 'temp4_std',
		'temp5_avg', 'temp5_max', 'temp5_min', 'temp5_std',
		'batt_a','batt_b'
	GM2:
		'table', 'year', 'doy', 'time',
		'temp1_avg', 'temp1_max', 'temp1_min', 'temp1_std',
		'temp2_avg', 'temp2_max', 'temp2_min', 'temp2_std',
		'temp3_avg', 'temp3_max', 'temp3_min', 'temp3_std',
		'temp4_avg', 'temp4_max', 'temp4_min', 'temp4_std',
		'temp5_avg', 'temp5_max', 'temp5_min', 'temp5_std',
		'batt_a','batt_b'



Resources:
	An example jupyter notebook for reading and plotting the data is available on github:
		https://github.com/spestana/snowex2020-snow-temp

Create a list of column headers according to the readme above (for “GM1” which we can read was the datalogger at snowpit 2S10)

column_headers = ['table', 'year', 'doy', 'time', # year, day of year, time of day (local time, UTC-7)
                  'rad_avg', 'rad_max', 'rad_min', 'rad_std', # radiometer surface temperature
                  'sb_avg', 'sb_max', 'sb_min', 'sb_std', # radiometer sensor body temperature (for calibration)
                  'temp1_avg', 'temp1_max', 'temp1_min', 'temp1_std', # temperature at 5 cm below snow surface
                  'temp2_avg', 'temp2_max', 'temp2_min', 'temp2_std', #               10 cm
                  'temp3_avg', 'temp3_max', 'temp3_min', 'temp3_std', #               15 cm
                  'temp4_avg', 'temp4_max', 'temp4_min', 'temp4_std', #               20 cm
                  'temp5_avg', 'temp5_max', 'temp5_min', 'temp5_std', #               30 cm
                  'batt_a','batt_b', # battery voltage data
                 ]

Open the file as a pandas data frame with read_csv

df = pd.read_csv(f'data/CR10X_GM1_final_storage_1.dat',
                 header = None, names = column_headers) 

# After the filepath we specify header=None because the file doesn't contain column headers, 
# then we specify names=column_headers to give our own names for each column.

We need to do some formatting of the data fields, but we can preview what we just loaded fist

df.head() # show the first 5 rows of the dataframe
Loading...

Data cleanup and formatting

# Create a zero-padded time string (e.g. for 9:30 AM we are changing '930' into '0930')
df['time_str'] = [('0' * (4 - len(str(df.time[i])))) + str(df.time[i]) for i in range(df.shape[0])]

# locate where rows have time_str == 2400 (midnight), and the whole column 'doy'
# where we are at midnight, we need to shift one day forward
df.loc[df['time_str'] == '2400','doy'] += 1

# and then change midnight from '2400' to '0000'
df.time_str.replace('2400', '0000', inplace=True)
/tmp/ipykernel_3843/748583221.py:9: FutureWarning: A value is trying to be set on a copy of a DataFrame or Series through chained assignment using an inplace method.
The behavior will change in pandas 3.0. This inplace method will never work because the intermediate object on which we are setting values always behaves as a copy.

For example, when doing 'df[col].method(value, inplace=True)', try using 'df.method({col: value}, inplace=True)' or df[col] = df[col].method(value) instead, to perform the operation inplace on the original object.


  df.time_str.replace('2400', '0000', inplace=True)

This function lets us convert year and day of year (the format that the datalogger uses) to a pandas datetime index:

def compose_date(years, months=1, days=1, weeks=None, hours=None, minutes=None,
                 seconds=None, milliseconds=None, microseconds=None, nanoseconds=None):
    '''Compose a datetime object from various datetime components. This clever solution is from:
        https://stackoverflow.com/questions/34258892/converting-year-and-day-of-year-into-datetime-index-in-pandas'''
    years = np.asarray(years) - 1970
    months = np.asarray(months) - 1
    days = np.asarray(days) - 1
    types = ('<M8[Y]', '<m8[M]', '<m8[D]', '<m8[W]', '<m8[h]',
             '<m8[m]', '<m8[s]', '<m8[ms]', '<m8[us]', '<m8[ns]')
    vals = (years, months, days, weeks, hours, minutes, seconds,
            milliseconds, microseconds, nanoseconds)
    return sum(np.asarray(v, dtype=t) for t, v in zip(types, vals)
               if v is not None)
# Create a datetime value from the date field and zero-padded time_str field, set this as our dataframe's index
df.index = compose_date(df['year'], 
                        days=df['doy'], 
                        hours=df['time_str'].str[:2],
                        minutes=df['time_str'].str[2:])

# Remove entries that are from table "102" (this contains datalogger battery information we're not interested in at the moment)
df = df[df.table != 102]

# drop the columns we no longer need
df.drop(columns=['table','year','doy','time','time_str','batt_a','batt_b'], inplace=True)

Inspect the contents

df.head()
Loading...

Make a simple plot of the data. We are interested in the variable rad_avg which is the average temperature measured by the radiometer over each 5 minute period.

plt.figure(figsize=(10,4))

# plot radiometer average temperature
df.rad_avg.plot(linestyle='-', marker='', markersize=1, c='k', label='Ground-based $T_s$')


# set axes limits
plt.ylim((-35,5))
plt.xlim((pd.Timestamp(2020,2,5,11,0),pd.Timestamp(2020,2,12,16,0)))

# add a legend to the plot
plt.legend()

# set axes labels
plt.ylabel(r'Temperature [$C\degree$]')
plt.xlabel('Time')

# add grid lines to the plot
plt.grid('on')

# set the plot title
plt.title('Snow Surface Temperature at Snow Pit 2S10');
<Figure size 1000x400 with 1 Axes>
Bonus plot: look at snow temperatures below the snow surface

Add the following to the above plot to add lines for temperature recorded at each depth interval below the snow surface:

# plot the snow temperature at each depth it was measured
df.temp1_avg.plot(linestyle='-', marker='.', markersize=1, c=[0.8,0.8,1], label='Ts @ -5 cm')
df.temp2_avg.plot(linestyle='-', marker='.', markersize=1, c=[0.6,0.6,1], label='Ts @ -10 cm')
df.temp3_avg.plot(linestyle='-', marker='.', markersize=1, c=[0.4,0.4,1], label='Ts @ -15 cm')
df.temp4_avg.plot(linestyle='-', marker='.', markersize=1, c=[0.2,0.2,1], label='Ts @ -20 cm')
df.temp5_avg.plot(linestyle='-', marker='.', markersize=1, c=[0,0,1], label='Ts @ -30 cm')

But then we want to focus on the date/time when our IR image was from, so zoom in on Feb 8th by changing our plot’s xlim (using pandas Timestamps for the x axis values).

plt.figure(figsize=(10,4))

# plot radiometer average temperature
df.rad_avg.plot(linestyle='-', marker='', markersize=1, c='k', label='Ground-based $T_s$')

# set axes limits
plt.ylim((-15,0)) # set some temperature y-axis limits for our plot
plt.xlim((pd.Timestamp(2020,2,8,6,0),pd.Timestamp(2020,2,8,20,0))) # zoom in to daytime hours on Feb. 8, 2020

# add a legend to the plot
plt.legend()

# set axes labels
plt.ylabel(r'Temperature [$C\degree$]')
plt.xlabel('Time')

# add grid lines to the plot
plt.grid('on')

# set the plot title
plt.title('Snow Surface Temperature at Snow Pit 2S10');
<Figure size 1000x400 with 1 Axes>

Compare Airborne IR against the “ground truth” snow surface temperature

What is the temperature at this point in the airborne IR image?

Use rioxarray’s clip function to extract the raster values that intersect with the point’s geometry. Because we have a point, this will return a single value for the pixel that overlaps this point.

# clip using our point's geometry
airborne_ir_point_temperature = airborne_ir.rio.clip(siteData_df.geometry)

# preview the result
airborne_ir_point_temperature
Loading...

Our result is a DataArray with a single data value at one set of x and y coordinates.

Add a 100 m radius buffer around this point and get the temperature from the airborne imagery for a larger area around the snow pit.

r = 100 # radius of the buffer in meters (this is in meters because we are working in a UTM coordinate reference system)

# create the buffered geometry
siteData_df_buffer = siteData_df.buffer(r)

# preview the resulting geometry, we should see this is a POLYGON now
siteData_df_buffer
0 POLYGON ((743176 4322689, 743175.518 4322679.1... 1 POLYGON ((743176 4322689, 743175.518 4322679.1... 2 POLYGON ((743176 4322689, 743175.518 4322679.1... 3 POLYGON ((743176 4322689, 743175.518 4322679.1... 4 POLYGON ((743176 4322689, 743175.518 4322679.1... dtype: geometry

What does this polygon look like when we plot it on top of the airborne IR image now?

fig, ax = plt.subplots(figsize=(10,10)) # create a new matplotlib figure, set the figure size
ax.set_aspect('equal') # set the aspect ratio to "equal"

# plot the airborne infrared image
airborne_ir.plot(cmap='magma', vmin=-10, vmax=10, ax=ax, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})

# plot the location of the snow pit of interest
siteData_df.plot(ax=ax, color='r', marker='x')

# plot the area of the buffer we made around the snow pit
siteData_df_buffer.plot(ax=ax, edgecolor='r', facecolor='none')

# set the same axes limits as above
ax.set_xlim((xmin-1000, xmax+1000)) # x axis limits to +/- 1 km from our point's "total bounds"
ax.set_ylim((ymin-1000, ymax+1000)) # y axis limits to +/- 1 km from our point's "total bounds"

# set axes labels
ax.set_xlabel('Eastings UTM 12N (m)')
ax.set_ylabel('Northings UTM 12N (m)')

# set plot title
ax.set_title('Location of Snow Pit 2S10, and 100 m radius buffer\nwith ASTER Band 14 TIR imagery');
<Figure size 1000x1000 with 2 Axes>

Clip the airborne IR raster again, now with our 200 m diameter polygon around the snow pit site.

# clip using our new geometry
airborne_ir_area_temperature = airborne_ir.rio.clip(siteData_df_buffer.geometry)

# preview the result
airborne_ir_area_temperature
Loading...

The result of clipping is again a DataArray, this time though it is 40x40. We can plot this to see what it looks like, and to see the distribution of temperatures in the area.

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10,5), tight_layout=True)

# plot the portion of the airborne TIR image we selected within the buffer area geometry
airborne_ir_area_temperature.plot(cmap='magma', vmin=-7, vmax=-4, ax=ax[0], 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})
ax[0].set_title('Airborne TIR image within\n100 m radius buffer (2S10)\n')
ax[0].set_aspect('equal')
ax[0].set_xlabel('Eastings UTM 12N (m)')
ax[0].set_ylabel('Northings UTM 12N (m)')
ax[0].set_xlim((xmin-150, xmax+150)) # x axis limits to +/- 150 m from our point's "total bounds"
ax[0].set_ylim((ymin-150, ymax+150)) # y axis limits to +/- 150 m from our point's "total bounds"

# plot the location of the snow pit of interest to the plot
siteData_df.plot(ax=ax[0], color='c', marker='x')

# plot the area of the buffer we made around the snow pit
siteData_df_buffer.plot(ax=ax[0], edgecolor='c', facecolor='none')

# plot a histogram of image temperature data within the buffer area geometry
airborne_ir_area_temperature.plot.hist(ax=ax[1],
                                       color='k', 
                                       zorder=1, # use zorder to make sure this plots below the point
                                       label='zonal $T_S$ histogram') 

# plot a vertical line for the single-pixel temperature we think is right at the snow pit
ax[1].axvline(airborne_ir_point_temperature, 
              color='c',linestyle='--',  # set color and style
              zorder=2, # use zorder to make sure this plots on top of the histogram
              label='$T_S$ single pixel') 

# plot a vertical line for the mean temperature within the buffer area geometry
ax[1].axvline(airborne_ir_area_temperature.mean(), 
              color='m',linestyle='--',  # set color and style
              zorder=2, # use zorder to make sure this plots on top of the histogram
              label='zonal mean $T_S$') 

ax[1].legend(loc='upper left') # add a legend
ax[1].set_xlim((-7,-4)) # set xlim to same values as colorbar in image plot
ax[1].set_ylim((0,400)) # set ylim
ax[1].set_title('Snow surface temperatures\nfrom Airborne TIR image (snow pit 2S10)')
ax[1].set_ylabel('Number of pixels');
<Figure size 1000x500 with 3 Axes>
What do these plots tell us about the surface temperature around the snow pit as measured by the airborne IR cameras?

  • Is the snow surface temperature more variable or uniform in this area?

    • Note that we changed the minimum and maximum values of our colorbar!

    • Try computing the standard deviation of temperatures in this area

  • How does our single pixel at the center compare with the rest of this area?

    • Try taking the mean or median of all temperatures in the area and compare against the single point. What is the difference?

  • See Lundquist et al., 2018 for an application of these methods with airborne and MODIS thermal infrared imagery

Plot the airborne IR temperature data on top of the ground-based timeseries

plt.figure(figsize=(10,4))

# plot radiometer average temperature
df.rad_avg.plot(linestyle='-', marker='', markersize=1, c='k', label='Ground-based $T_s$')

# plot the mean airborne IR temperature from the area around the snow pit:
plt.plot(airborne_ir_timestamp, airborne_ir_area_temperature.mean(),
         marker='o', c='r', linestyle='none',
         label='Airborne IR mean $T_s$ for 100 m radius area')

# plot an error bar showing the maximum and minimum airborne IR temperature around the snow pit
plt.errorbar(airborne_ir_timestamp, airborne_ir_area_temperature.mean(),
             yerr=[[airborne_ir_area_temperature.mean()-airborne_ir_area_temperature.min()], 
                   [airborne_ir_area_temperature.max()-airborne_ir_area_temperature.mean()]],
            capsize=3, fmt='none', ecolor='r',
            label='Airborne IR $T_s$ range for 100 m radius area')


# set axes limits
plt.ylim((-15,0))
plt.xlim((pd.Timestamp(2020,2,8,6,0),pd.Timestamp(2020,2,8,20,0))) # zoom in to daytime hours on Feb. 8, 2020

# add a legend to the plot
plt.legend()

# set axes labels
plt.ylabel(r'Temperature [$C\degree$]')
plt.xlabel('Time')

# add grid lines to the plot
plt.grid('on')

# set the plot title
plt.title('Snow Surface Temperature at Snow Pit 2S10');
/tmp/ipykernel_3843/1336368688.py:7: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  plt.plot(airborne_ir_timestamp, airborne_ir_area_temperature.mean(),
<Figure size 1000x400 with 1 Axes>
Continuing with this analysis:

In the above plot we’ve added “error bars” to represent the full range of temperatures within the 100 m radius area around the snow pit.

  • Is this a fair comparison?

  • Should we make the area smaller based on our confidence in the image’s geolocation accuracy (10-15 m)?

  • What is the difference between the “ground truth” data and the airborne IR data at this point in time?

Part 2: Satellite IR remote sensing obsevations

Satellite IR imagery with ASTER

Advantages of satellite IR images: We don’t always have airplanes with IR cameras flying around. Satellites can provide images at more regular intervals for long-term studies, and can see areas that are difficult to access on the ground or by air.

What might be some diadvantages of satellite IR imagery compared to airborne IR imagery?

  • Lower spatial resolution because they’re further away from the Earth’s surface

  • Mixed pixel problem: with lower image resolutions, each pixel contains a more heterogeneous mixtures of surfaces and temperatures, meaning that temperature information is more blurred together

  • What other disadvantages can you think of?

  • How might these differences (advantages or disadvantages) change the type of research questions you can investigate?

For this tutorial, we will look at an image from NASA’s Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imager, which is onboard the Terra satellite along with a MODIS imager. We can compare an ASTER IR image of Grand Mesa that was taken at roughly the same time as the airborne IR image. The ASTER image we will be working with is from ASTER’s band 14 which is sensitive to radiance in the 10.95-11.65 µm wavelength range.

Load an ASTER geotiff that we’ve downloaded for this tutorial, and inspect its contents.

aster_ir = rioxarray.open_rasterio(S3_BASE_URL + 'AST_L1T_00302082020180748_20200209065849_17218_ImageData14.tif')

Inspect the ASTER file we just opened

aster_ir
Loading...

What is its CRS?

aster_ir.rio.crs
CRS.from_wkt('PROJCS["WGS 84 / UTM zone 12N",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",-111],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","32612"]]')
# Reproject this ASTER image into our common coordinate system
aster_ir = aster_ir.rio.reproject('EPSG:26912') # overwrite itself with new reprojected data array
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/AST_L1T_00302082020180748_20200209065849_17218_ImageData14.tif.msk: 403
Warning 1: HTTP response code on https://snowex-tutorials.s3.us-west-2.amazonaws.com/thermal-ir/AST_L1T_00302082020180748_20200209065849_17218_ImageData14.tif.MSK: 403

When was this image taken?

aster_ir_timestamp = pd.Timestamp(2020,2,8,18,7,48) - pd.Timedelta(hours=7)

Plot the image. What are the units of the values on the colorbar?

aster_ir.plot()
<Figure size 640x480 with 2 Axes>

It’s necessary to read the product documentation to understand what we are looking at here.

We are using the ASTER Level 1 Precision Terrain Corrected Registered At-Sensor Radiance (AST_L1T) product. Product documentation is available here. Also helpful is the ** (AST_L1B) product documentation here from which AST_L1T is derived.

The values here are stored as scaled “digital number” (DN) values rather than the actual radiance values. The product documentation also provides information about how to unscale these values back into radiance units, and from radiance to brightness temperature, using laboratory calibrated constants.

I’ve written two functions here to do this unit conversion for the five ASTER TIR bands in two steps (DN to radiance, radiance to brightness temperature). The function takes as its arguments the DN or radiance values respectively, and the band number (in our case band number 14, not the band wavelengths).

def tir_dn2rad(DN, band):
    '''Convert AST_L1T Digital Number values to At-Sensor Radiance for the TIR bands (bands 10-14).'''
    ucc = [6.822e-3, 6.780e-3, 6.590e-3, 5.693e-3, 5.225e-3]
    rad = (DN-1.) * ucc[band-10]
    return rad

def tir_rad2tb(rad, band):
    '''Convert AST_L1T At-Sensor Radiance to Brightness Temperature [K] for the TIR bands (bands 10-14).'''
    k1 = [3047.47, 2480.93, 1930.80, 865.65, 649.60]
    k2 = [1736.18, 1666.21, 1584.72,1349.82, 1274.49]
    tb = k2[band-10] /  np.log((k1[band-10]/rad) + 1)
    return tb

Use the above functions to convert from DN to radiance, radiance to brightness temperature (assume and emissivity of 1 for all surfaces, note that the airborne imagery also assumed emissivity of 1).

Then convert from degeees K to degrees C by subtracting 273.15.

aster_band14_rad = tir_dn2rad( aster_ir, band=14 ) # convert from DN to radiance
aster_band14_tb_k = tir_rad2tb( aster_band14_rad, band=14 ) # convert from radiance to brightness temperature (K)
aster_band14_tb_c = aster_band14_tb_k - 273.15 # convert from K to C

# Note that an "invalid value encountered..." warning may pop up here. This is because the above function tries to take the log of "nan" values that are outside the imaged area
# we can ignore this warning and proceed
/home/runner/micromamba/envs/snow-observations-cookbook-test/lib/python3.13/site-packages/xarray/computation/apply_ufunc.py:818: RuntimeWarning: invalid value encountered in log
  result_data = func(*input_data)

During this unit conversion, xarray dropped the coordinate reference system attributes, so here we add the original crs to the new ASTER degrees celsius dataarray.

aster_band14_tb_c.rio.set_crs(aster_ir.rio.crs, inplace=True);
/tmp/ipykernel_3843/3871339668.py:1: FutureWarning: It is recommended to use 'rio.write_crs()' instead. 'rio.set_crs()' will likelybe removed in a future release.
  aster_band14_tb_c.rio.set_crs(aster_ir.rio.crs, inplace=True);

Plot our image again, this time setting our colorscale and colorbar values. We should see “realistic” surface temperature values in degrees C now.

fig, ax = plt.subplots(figsize=(15,10))
ax.set_aspect('equal')
aster_band14_tb_c.plot(ax=ax,
                       cmap='magma', 
                       vmin=-20, vmax=20, # note that we have a wider temperature range on our colorbar for this image, -20 to +20 C instead of -10 to +10 C
                       cbar_kwargs={'label': r'Temperature $\degree C$'})
<Figure size 1500x1000 with 2 Axes>

Plot ASTER next to Airborne IR to see spatial resolution differences

fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(12,6), tight_layout=True)

aster_band14_tb_c.plot(ax=axs[0], 
                       cmap='magma', 
                       vmin=-10, vmax=10, 
                       cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[0].set_title('ASTER IR')
axs[0].set_aspect('equal')


airborne_ir.plot(ax=axs[1], 
                 cmap='magma', 
                 vmin=-10, vmax=10, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[1].set_title('Airborne IR')
axs[1].set_aspect('equal')

# for each subplot, do the following:
for ax in axs: 
    ax.set_aspect('equal') # set the aspect ratio to "equal"
    
    # plot the location of the snow pit of interest
    siteData_df.plot(ax=ax, color='r', marker='x')
    # plot the area of the buffer we made around the snow pit
    siteData_df_buffer.plot(ax=ax, edgecolor='r', facecolor='none')
    
    # set axes labels and limits
    ax.set_xlabel('Eastings UTM 12N (m)')
    ax.set_xlim((735000, 760000))
    ax.set_ylabel('Northings UTM 12N (m)')
    ax.set_ylim((4320000, 4325000))
<Figure size 1200x600 with 4 Axes>

Make another plot to zoom in on snow pit 2S10:

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12,4), tight_layout=True)

aster_band14_tb_c.plot(ax=axs[0], 
                       cmap='magma', 
                       vmin=-10, vmax=10, 
                       cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[0].set_title('ASTER IR')
axs[0].set_aspect('equal')


airborne_ir.plot(ax=axs[1], 
                 cmap='magma', 
                 vmin=-10, vmax=10, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[1].set_title('Airborne IR')
axs[1].set_aspect('equal')

# for each subplot, do the following:
for ax in axs: 
    ax.set_aspect('equal') # set the aspect ratio to "equal"
    
    # plot the location of the snow pit of interest
    siteData_df.plot(ax=ax, color='r', marker='x')
    # plot the area of the buffer we made around the snow pit
    siteData_df_buffer.plot(ax=ax, edgecolor='r', facecolor='none')
    
    # set axes labels
    ax.set_xlabel('Eastings UTM 12N (m)')
    ax.set_ylabel('Northings UTM 12N (m)')
    # set the same axes limits as above
    ax.set_xlim((xmin-1000, xmax+1000)) # x axis limits to +/- 1 km from our point's "total bounds"
    ax.set_ylim((ymin-1000, ymax+1000)) # y axis limits to +/- 1 km from our point's "total bounds"
<Figure size 1200x400 with 4 Axes>

Get the temperature of the ASTER pixel at the snow pit point using rioxarray clip.

# First clip to the single point
aster_band14_tb_c_point_temperature = aster_band14_tb_c.rio.clip(siteData_df.geometry)

# Second clip to the 100 m radius buffered area
aster_band14_tb_c_area_temperature = aster_band14_tb_c.rio.clip(siteData_df_buffer.geometry)

# preview the result
aster_band14_tb_c_area_temperature
Loading...

How many ASTER pixels in the area did we select?

Plot the clipped area:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10,5), tight_layout=True)

# plot the portion of the airborne TIR image we selected within the buffer area geometry
aster_band14_tb_c_area_temperature.plot(cmap='magma', vmin=-7, vmax=-4, ax=ax[0], 
                 cbar_kwargs={'label': 'Temperature $\degree C$'})
ax[0].set_title('ASTER TIR image within\n100 m radius buffer (2S10)\n')
ax[0].set_aspect('equal')
ax[0].set_xlabel('Eastings UTM 12N (m)')
ax[0].set_ylabel('Northings UTM 12N (m)')
ax[0].set_xlim((xmin-150, xmax+150)) # x axis limits to +/- 150 m from our point's "total bounds"
ax[0].set_ylim((ymin-150, ymax+150)) # y axis limits to +/- 150 m from our point's "total bounds"

# plot the location of the snow pit of interest to the plot
siteData_df.plot(ax=ax[0], color='r', marker='x')

# plot the area of the buffer we made around the snow pit
siteData_df_buffer.plot(ax=ax[0], edgecolor='r', facecolor='none')

# plot a histogram of image temperature data within the buffer area geometry
aster_band14_tb_c_area_temperature.plot.hist(ax=ax[1], color='k');
ax[1].set_xlim((-7,-4)) # set xlim to same values as colorbar in image plot
ax[1].set_title('Histogram of temperatures from ASTER TIR image\n100 m radius buffer (2S10)')
ax[1].set_ylabel('Number of pixels');
<>:5: SyntaxWarning: invalid escape sequence '\d'
<>:5: SyntaxWarning: invalid escape sequence '\d'
/tmp/ipykernel_3843/3974720511.py:5: SyntaxWarning: invalid escape sequence '\d'
  cbar_kwargs={'label': 'Temperature $\degree C$'})
<Figure size 1000x500 with 3 Axes>

Plot the ASTER IR temperature data on top of the ground-based timeseries and airborne IR temperature

plt.figure(figsize=(10,4))

# plot radiometer average temperature
df.rad_avg.plot(linestyle='-', marker='', markersize=1, c='k', label='Ground-based $T_s$')

# plot the mean airborne IR temperature from the area around the snow pit:
plt.plot(airborne_ir_timestamp, airborne_ir_area_temperature.mean(),
         marker='o', c='r', linestyle='none',
         label='Airborne IR mean $T_s$ for 100 m radius area')
# plot an error bar showing the maximum and minimum airborne IR temperature around the snow pit
plt.errorbar(airborne_ir_timestamp, airborne_ir_area_temperature.mean(),
             yerr=[[airborne_ir_area_temperature.mean()-airborne_ir_area_temperature.min()], 
                   [airborne_ir_area_temperature.max()-airborne_ir_area_temperature.mean()]],
            capsize=3, fmt='none', ecolor='r',
            )#label='Airborne IR $T_s$ range for 100 m radius area')

# plot the mean ASTER IR temperature from the area around the snow pit:
plt.plot(aster_ir_timestamp, aster_band14_tb_c_area_temperature.mean(),
         marker='o', c='b', linestyle='none',
        label='ASTER IR mean $T_s$ for 100 m radius area')
# plot an error bar showing the maximum and minimum ASTER IR temperature around the snow pit
plt.errorbar(aster_ir_timestamp, aster_band14_tb_c_area_temperature.mean(),
             yerr=[[aster_band14_tb_c_area_temperature.mean()-aster_band14_tb_c_area_temperature.min()], 
                   [aster_band14_tb_c_area_temperature.max()-aster_band14_tb_c_area_temperature.mean()]],
            capsize=3, fmt='none', ecolor='b',
            )#label='ASTER IR $T_s$ range for 100 m radius area')


# set axes limits
plt.ylim((-15,0))
plt.xlim((pd.Timestamp(2020,2,8,6,0),pd.Timestamp(2020,2,8,20,0))) # zoom in to daytime hours on Feb. 8, 2020

# add a legend to the plot
plt.legend()

# set axes labels
plt.ylabel('Temperature [$C\degree$]')
plt.xlabel('Time')

# add grid lines to the plot
plt.grid('on')

# set the plot title
plt.title('Snow Surface Temperature at Snow Pit 2S10');
<>:37: SyntaxWarning: invalid escape sequence '\d'
<>:37: SyntaxWarning: invalid escape sequence '\d'
/tmp/ipykernel_3843/4068806060.py:37: SyntaxWarning: invalid escape sequence '\d'
  plt.ylabel('Temperature [$C\degree$]')
/tmp/ipykernel_3843/4068806060.py:7: UserWarning: This axis already has a converter set and is updating to a potentially incompatible converter
  plt.plot(airborne_ir_timestamp, airborne_ir_area_temperature.mean(),
<Figure size 1000x400 with 1 Axes>

Bonus activity: comparing two thermal IR rasters

How does the finer spatial resolution airborne IR image compare with the coarser resolution ASTER IR image?

One way to compare these two images is to “upscale” the finer resolution airborne IR image to the same spatial resolution as ASTER.

We can use the rioxarray reproject_match function to do this. (also see this example)

Note that this function has multiple options for how we want to resample the image data to the new spatial resolution. We will use resampling=5 which corresponds to taking the mean value.

airborne_ir_repr = airborne_ir.rio.reproject_match(aster_ir, resampling=5)

Preview the result.

  • What are the image dimensions of the reprojected airborne IR image compared with the original?

airborne_ir_repr
Loading...

Plot the ASTER image, original airborne IR image, and resampled airborne IR image next to each other to visualize these spatial resolution differences.

fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(15,4), tight_layout=True)

# ASTER IR image
aster_band14_tb_c.plot(ax=axs[0], 
                       cmap='magma', 
                       vmin=-10, vmax=10, 
                       cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[0].set_title('ASTER IR')
axs[0].set_aspect('equal')

# Original airborne IR image
airborne_ir.plot(ax=axs[1], 
                 cmap='magma', 
                 vmin=-10, vmax=10, 
                 cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[1].set_title('Airborne IR')
axs[1].set_aspect('equal')

# Resampled airborne IR image
airborne_ir_repr.plot(ax=axs[2], 
                      cmap='magma', 
                      vmin=-10, vmax=10, 
                      cbar_kwargs={'label': r'Temperature $\degree C$'})
axs[2].set_title('Airborne IR resampled (mean)\nto ASTER resolution (90m)')
axs[2].set_aspect('equal')


# for each subplot, do the following:
for ax in axs: 
    ax.set_aspect('equal') # set the aspect ratio to "equal"
    
    # plot the location of the snow pit of interest
    siteData_df.plot(ax=ax, color='r', marker='x')
    # plot the area of the buffer we made around the snow pit
    siteData_df_buffer.plot(ax=ax, edgecolor='r', facecolor='none')
    
    # set axes labels
    ax.set_xlabel('Eastings UTM 12N (m)')
    ax.set_ylabel('Northings UTM 12N (m)')
    # set the same axes limits as above
    ax.set_xlim((xmin-1000, xmax+1000)) # x axis limits to +/- 1 km from our point's "total bounds"
    ax.set_ylim((ymin-1000, ymax+1000)) # y axis limits to +/- 1 km from our point's "total bounds"
<Figure size 1500x400 with 6 Axes>

Finally, compute the difference between the ASTER IR image and resampled airborne IR image, then plot the result:

# Subtract reprojected airborne IR image from ASTER IR image
difference_image = aster_band14_tb_c - airborne_ir_repr
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6,5), tight_layout=True)

# Difference image
difference_image.plot(ax=ax, 
                       cmap='RdBu_r', 
                       vmin=-10, vmax=10, 
                       cbar_kwargs={'label': 'Temperature $\degree C$'})
ax.set_title('Temperature Difference\n(ASTER IR - Reprojected Airborne IR)\n')
ax.set_aspect('equal')

# plot the location of the snow pit of interest
siteData_df.plot(ax=ax, color='r', marker='x')
# plot the area of the buffer we made around the snow pit
siteData_df_buffer.plot(ax=ax, edgecolor='r', facecolor='none')

# set axes labels
ax.set_xlabel('Eastings UTM 12N (m)')
ax.set_ylabel('Northings UTM 12N (m)')
# set the same axes limits as above
ax.set_xlim((xmin-1000, xmax+1000)) # x axis limits to +/- 1 km from our point's "total bounds"
ax.set_ylim((ymin-1000, ymax+1000)); # y axis limits to +/- 1 km from our point's "total bounds"
<>:7: SyntaxWarning: invalid escape sequence '\d'
<>:7: SyntaxWarning: invalid escape sequence '\d'
/tmp/ipykernel_3843/1907299178.py:7: SyntaxWarning: invalid escape sequence '\d'
  cbar_kwargs={'label': 'Temperature $\degree C$'})
<Figure size 600x500 with 2 Axes>

Next steps:

Project ideas with these data:

  • Compare more snow pit temperature data (using snowexsql queries) against airborne and satellite IR imagery

  • Investigate spatial patterns of snow temperature from open areas to forested areas on the mesa, differences between ASTER and airborne IR imagery

  • Improved data visualization using hvplot or something similar to create interactive plots of IR and/or visible imagery

  • Compare thermal infrared and SAR imagery, or snow temperature observations and snow model outputs

Data access/download and pre-processing

See the jupyter notebook thermal-ir-data-download.ipynb for more details about data access methods through the NASA EarthData API, and pre-processing ASTER data into geotiff images.

Additional learning resources:

Open an airborne TIR mosaic NetCDF file

# Open airborne TIR mosaic NetCDF file
with fsspec.open(S3_BASE_URL + 'SNOWEX2020_IR_PLANE_2020Feb08_mosaicked_APLUW.nc') as f:
    ds = xr.open_dataset(f, decode_times=False)
    ds = ds.load()

Inspect the dataset and its dimensions

# Preview the dataset
ds
Loading...
# Take a look at the dimensions in the dataset
ds.dims
FrozenMappingWarningOnValuesAccess({'pass': 17, 'easting, x': 4171, 'northing, y': 3581, 'na': 1})

There is an extra dimension (“na”) that we will want to drop, we will want to rename some of the dims, assign coordinates to those dims, and add a coordinate reference system for plotting.

# Drop the extra "na" dimension from E_UTM, N_UTM, and time
ds['E_UTM'] = ds['E_UTM'].isel(na=0, drop=True)
ds['N_UTM'] = ds['N_UTM'].isel(na=0, drop=True)
ds['time'] = ds['time'].isel(na=0, drop=True)

Rename the dimensions to some easier to use names

# Rename dims
ds = ds.rename({"pass" : "time", 
                "easting, x" : "easting", 
                "northing, y" : "northing"})

This NetCDF file was generated in MATLAB, and the dates/times are in an epoch format. Use utcfromtimestamp() and isoformat() to convert and reformat into a more convenient format.

# Decode matlab (epoch) format times
utctime = [datetime.utcfromtimestamp(this_time).isoformat() for this_time in ds.time.values]
/tmp/ipykernel_3843/3153572546.py:2: DeprecationWarning: datetime.datetime.utcfromtimestamp() is deprecated and scheduled for removal in a future version. Use timezone-aware objects to represent datetimes in UTC: datetime.datetime.fromtimestamp(timestamp, datetime.UTC).
  utctime = [datetime.utcfromtimestamp(this_time).isoformat() for this_time in ds.time.values]

Assign and then transpose coordinates in our dataset

# Assign coordinates to the "northing",  "easting", and "time" dimensions
ds = ds.assign_coords({"time": utctime, "northing": ds.N_UTM, "easting": ds.E_UTM})
# Transpose coords
ds = ds.transpose("time", "northing", "easting")

Set spatial dimensions then define which coordinate reference system the spatial dimensions are in

# Write the coordinate reference system for the spatial dims with rioxarray
# https://github.com/corteva/rioxarray/issues/379
ds.rio.write_crs('epsg:32612', inplace=True,
                ).rio.set_spatial_dims('easting', 'northing', inplace=True,
                                      ).rio.write_coordinate_system(inplace=True)
Loading...
References
  1. Lundquist, J. D., Chickadel, C., Cristea, N., Currier, W. R., Henn, B., Keenan, E., & Dozier, J. (2018). Separating snow and forest temperatures with thermal infrared remote sensing. Remote Sensing of Environment, 209, 764–779. 10.1016/j.rse.2018.03.001
Observations
GPR and Lidar
Observations
Time-lapse Cameras