<jupyter_start><jupyter_text>Downloading datasets for the thermal IR tutorialThis is an optional notebook to demonstrate a method for accessing SnowEx data through the NASA EarthData API. Alternatively (as seen within the [thermal-ir-tutorial.ipynb](thermal-ir-tutorial.ipynb) notebook) you can access all data needed for the tutorial by running the following command in a terminal:`aws s3 sync --quiet s3://snowex-data/tutorial-data/thermal-ir/ /tmp/thermal-ir/`:::{admonition} Learning Objectives**At the conclusion of this tutorial, you will be able to:**- download datasets through the NASA Earthdata API- use data conversion tools to create geotiff files:::<jupyter_code># Import the packages we'll need
import pandas as pd
from earthdata_api import earthdata_granule_search # for interfacing with the NASA Earthdata API<jupyter_output><empty_output><jupyter_text>We are going to use a custom function (`earthdata_api.earthdata_granule_search()`) which can return a list of download URLs for files we search for. We can then use a utility like wget (with our EarthData credentials) to download the files from that list of URLs. Alternatively, we could search on the [EarthData](https://earthdata.nasa.gov/) or [NSIDC](https://nsidc.org/) website itself and download through the web interface. **First, set up some of the search parameters we'll use for all the datasets we want to find:*** a bounding box for the Grand Mesa area, * and a date/time that we want to look for. When defining a bounding box for the study area we want to look at, Earthdata will return anything that overlaps this spot.What we are specifying  below is a a tiny area, basicallcy a point search, but we can expand to a larger area if needed.Format the bounding box as a list of coordinates like `bounding_box = [ lower left longitude, lower left latitude, upper right longitude, upper right latitude ]`<jupyter_code>bounding_box = [-108.19, 39.03, -108.18, 39.04] # a little box on top of Grand Mesa<jupyter_output><empty_output><jupyter_text>We are going to look at one specific day, but we could search multiple dates/times. Those dates/times don't need to be a continuous time period, we could instead provide a list of dates/times. (That was my original motivation to write this function). Use either [pandas.Timestamp](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Timestamp.html), [datetime.datetime](https://docs.python.org/3/library/datetime.htmldatetime.datetime), or [numpy.datetime64](https://numpy.org/doc/stable/reference/arrays.datetime.html) objects.<jupyter_code>timestampsUTC = pd.Timestamp('2020-2-08') # February 8, 2020<jupyter_output><empty_output><jupyter_text>---**Now we can find the first dataset we want, [ASTER](https://asterweb.jpl.nasa.gov/) satellite images.** Specificallfy we are going to download the AST_L1T product. Looking its [product documentation](https://lpdaac.usgs.gov/products/ast_l1tv003/) we see that it is available either as an HDF file with all bands, or as tif files with select visible or thermal bands.We want to download the hdf file, which we'll later convert to separate geotiffs for each band.<jupyter_code># Information for the ASTER product we want
product="AST_L1T" # product name
version="003" # product version
ext="hdf" # geotiff files only<jupyter_output><empty_output><jupyter_text>To make sure we find the images from this day, we need to specify a window in time to search, which starts at our `timestampsUTC` time we already defined.To define this time window we specify two things, the size of the time window from our start time, and units of that size value.<jupyter_code>search_time_window=1 # specify that we want to have a time window of 1
search_time_window_units="d" # specify time window units of Days, which together gives us 1 day<jupyter_output><empty_output><jupyter_text>We also want to save the list of data access URLs so we can use wget to download them<jupyter_code>list_filepath="/tmp/thermal-ir/aster_download_list.txt"
# make sure this directory exists
import os
if not os.path.exists("/tmp/thermal-ir"):
    os.mkdir("/tmp/thermal-ir")<jupyter_output><empty_output><jupyter_text>**Finally, make the API call to get the list of files that match our search**<jupyter_code>url_list = earthdata_granule_search(product, version, 
                                    bounding_box, timestampsUTC, 
                                    search_time_window, search_time_window_units, 
                                    list_filepath, ext)
# we can print the list of files, but these are also saved to a file
print(url_list)<jupyter_output><empty_output><jupyter_text>**Now we can open a terminal and use wget to download the files that are listed.**Open a terminal and navigate to the data directory for this tutorial`cd /tmp/thermal-ir/`Use wget to download the files listed in aster_download_list.txt. This requires that you put your Earthdata username wherer it says YOUR_USERNAME. This will also prompt you for your Earthdata password.`wget --http-user=YOUR_USERNAME --ask-password --keep-session-cookies --auth-no-challenge=on -c -i aster_download_list.txt` **Convert hdf to geotiffs with ASTERL1T_hdf2tir.py script.** When we run this utility, the output data will be at-sensor radiance stored as Digital Numbers (DN) as uint16 data types.The ASTERL1T_hdf2tir.py is from:*NASA LP DAAC. (2020). Reformat and Georeference ASTER L1T HDF Files Data Prep Script. Ver. 1. NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, USA. Accessed July 27, 2020. [https://git.earthdata.nasa.gov/projects/LPDUR/repos/aster-l1t/browse](https://git.earthdata.nasa.gov/projects/LPDUR/repos/aster-l1t/browse)*In the terminal, navigate to the directory containing the ASTERL1T_hdf2tir.py script`cd ~/website/book/tutorials/thermal-ir/ast-l1t/`Then run the script as follows to generate geotiff files for each ASTER band`python ASTERL1T_hdf2tif.py /tmp/thermal-ir/`You should now see the AST_L1T geotiff files in the data directory! ---**Download the second dataset we want, [ground-based temperature measurements](https://nsidc.org/data/SNEX20_VPTS_Raw/versions/1)**<jupyter_code>product = "SNEX20_VPTS_Raw"
version = "1"
list_filepath="/tmp/thermal-ir/snow_temperature_download_list.txt"
url_list = earthdata_granule_search(product, version, 
                                    bounding_box, timestampsUTC, 
                                    search_time_window, search_time_window_units, 
                                    list_filepath)
print(url_list)<jupyter_output><empty_output><jupyter_text>**Use wget to download these data**`wget --http-user=YOUR_USERNAME --ask-password --keep-session-cookies --auth-no-challenge=on -c -i snow_temperature_download_list.txt` Then unzip the folder we just downloaded (`-q` for "quiet mode" so it doesn't list every file as it unzips them, `-d` for the directory we want to unzip the files into)`unzip -q -d /tmp/thermal-ir/SNEX20_VPTS_Raw SNEX20_VPTS_Raw.zip`And you should see the folder "Level-0" with the raw ground-based snow temperature data inside now in your data directoryThe file we want is `/tmp/thermal-ir/SNEX20_VPTS_Raw/Level-0/snow-temperature-timeseries/CR10X_GM1_final_storage_1.dat`, this is a comma-delimited text file from the datalogger at snow pit 2S10--- ```{warning} Note:The airborne TIR imagery from SnowEx 2020 is not yet publicly available through NSIDC, we will work with a sample image provided in this tutorial.```An alternate method for re-formatting the airborne IR mosaic NetCDF files is below. (You must first download the sample file with `aws s3 sync --quiet s3://snowex-data/tutorial-data/thermal-ir/ /tmp/thermal-ir/` in the terminal)<jupyter_code>!aws s3 sync --quiet s3://snowex-data/tutorial-data/thermal-ir/ /tmp/thermal-ir/
import xarray as xr
import rioxarray
import datetime<jupyter_output><empty_output><jupyter_text>Open an airborne TIR mosaic NetCDF file<jupyter_code># Open airborne TIR mosaic NetCDF file
ds = xr.open_dataset('/tmp/thermal-ir/SNOWEX2020_IR_PLANE_2020Feb08_mosaicked_APLUW.nc',
                     decode_times=False) # set decode_times to False<jupyter_output><empty_output><jupyter_text>Inspect the dataset and its dimensions<jupyter_code># Preview the dataset
ds
# Take a look at the dimensions in the dataset
ds.dims<jupyter_output><empty_output><jupyter_text>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.<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>Rename the dimensions to some easier to use names<jupyter_code># Rename dims
ds = ds.rename({"pass" : "time", 
                "easting, x" : "easting", 
                "northing, y" : "northing"})<jupyter_output><empty_output><jupyter_text>This NetCDF file was generated in MATLAB, and the dates/times are in an epoch format. Use [utcfromtimestamp()](https://docs.python.org/3/library/datetime.htmldatetime.datetime.utcfromtimestamp) and [isoformat()](https://docs.python.org/3/library/datetime.htmldatetime.date.isoformat) to convert and reformat into a more convenient format.<jupyter_code># Decode matlab (epoch) format times
utctime = [datetime.datetime.utcfromtimestamp(this_time).isoformat() for this_time in ds.time.values]<jupyter_output><empty_output><jupyter_text>Assign and then transpose coordinates in our dataset<jupyter_code># 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")<jupyter_output><empty_output><jupyter_text>Set spatial dimensions then define which coordinate reference system the spatial dimensions are in<jupyter_code># Set spatial dimensions with rioxarray
ds.rio.set_spatial_dims('easting', 'northing', inplace=True);
# Write the coordinate reference system for the spatial dims with rioxarray
ds.rio.write_crs('epsg:32612', inplace=True);<jupyter_output><empty_output><jupyter_text>Inspect the final, re-formatted dataset<jupyter_code>ds<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09024925986056728, 'fraction_numerical': 0.01165122719893038, 'mean_word_length': 4.923076923076923, 'pattern_counts': {'":': 3, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 12, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '40018', 'license': 'MIT', 'n_tokens_mistral': 3337, 'n_tokens_neox': 3140, 'n_words': 1112, 'path': 'book/tutorials/thermal-ir/thermal-ir-data-download.ipynb', 'repo_name': 'meganmason/website', 'size': '10471'}