{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Working with ECOSTRESS Evapotranspiration Data\n",
"### This tutorial demonstrates how to work with the ECOSTRESS Evapotranspiration PT-JPL Daily L3 Global 70m Version 1 ([ECO3ETPTJPL.001](https://doi.org/10.5067/ECOSTRESS/ECO3ETPTJPL.001)) data product in Python.\n",
"The Land Processes Distributed Active Archive Center (LP DAAC) distributes the Ecosystem Spaceborne Thermal Radiometer Experiment on Space Station ([ECOSTRESS](https://ecostress.jpl.nasa.gov/)) data products. The ECOSTRESS mission is tasked with measuring the temperature of plants to better understand how much water plants need and how they respond to stress. ECOSTRESS products are archived and distributed in the HDF5 file format as swath-based products.\n",
"\n",
"In this tutorial, you will use Python to perform a swath to grid conversion to project the swath data on to a grid with a defined coordinate reference system (CRS), compare ECOSTRESS data with ground-based AmeriFlux flux tower observations, and export science dataset (SDS) layers as GeoTIFF files that can be loaded into a GIS and/or Remote Sensing software program. \n",
"***\n",
"### Example: Converting a swath ECO3ETPTJPL.001 HDF5 file into a GeoTIFF with a defined CRS and comparing ECOSTRESS Evapotranspiration (ET) with ground-based ET observations from an AmeriFlux flux tower location in California. \n",
"#### Data Used in the Example: \n",
"- Data Product: ECOSTRESS Evapotranspiration PT-JPL Daily L3 Global 70m Version 1 ([ECO3ETPTJPL.001](https://doi.org/10.5067/ECOSTRESS/ECO3ETPTJPL.001)) \n",
" - Science Dataset (SDS) layers: \n",
" - ETinst \n",
" - ETinstUncertainty \n",
"- Data Product: ECOSTRESS Geolocation Daily L1B Global 70m Version 1 ([ECO1BGEO.001](https://doi.org/10.5067/ECOSTRESS/ECO1BGEO.001)) \n",
" - Science Dataset (SDS) layers: \n",
" - Latitude \n",
" - Longitude \n",
"- Data Product: AmeriFlux Ground Observations for [Flux Tower US-CZ3](https://ameriflux.lbl.gov/sites/siteinfo/US-CZ3): Sierra Critical Zone, Sierra Transect, Sierran Mixed Conifer, P301\n",
" - Variables: \n",
" - Latent Heat (W/m$^{2}$) \n",
"*** \n",
"# Topics Covered:\n",
"1. **Getting Started** \n",
" 1a. Import Packages \n",
" 1b. Set Up the Working Environment \n",
" 1c. Retrieve Files \n",
"2. **Importing and Interpreting Data** \n",
" 2a. Open an ECOSTRESS HDF5 File and Read File Metadata \n",
" 2b. Subset SDS Layers \n",
"3. **Performing Swath2grid Conversion** \n",
" 3a. Import Geolocation File \n",
" 3b. Define Projection and Output Grid \n",
" 3c. Read SDS Metadata \n",
" 3d. Perform K-D Tree Resampling \n",
" 3e. Basic Image Processing \n",
"4. **Exporting Results** \n",
" 4a. Set Up a Dictionary \n",
" 4b. Define CRS and Export as GeoTIFFs \n",
"5. **Combining ECOSTRESS and AmeriFlux Tower Data** \n",
" 5a. Loading Tables with Pandas \n",
" 5b. Locate ECOSTRESS Pixel from Lat/Lon Coordinates \n",
"6. **Visualizing Data** \n",
" 6a. Create Colormap \n",
" 6b. Plot ET Data \n",
" 6c. Exporting an Image \n",
"7. **Comparing Observations** \n",
" 7a. Calculate Distribution of ECOSTRESS Data \n",
" 7b. Visualize Ground Observations \n",
" 7c. Combine ECOSTRESS and Ground Observations \n",
"\n",
"***\n",
"# Before Starting this Tutorial: \n",
"#### If you are simply looking to batch process/perform the swath2grid conversion for ECOSTRESS files, be sure to check out the [ECOSTRESS Swath to Grid Conversion Script](https://git.earthdata.nasa.gov/projects/LPDUR/repos/ecostress_swath2grid/browse).\n",
"\n",
"NOTE: This tutorial was developed specifically for the ECOSTRESS Evapotranspiration PT-JPL Level 3, Version 1 HDF5 files and will need to be adapted to work with other ECOSTRESS products. \n",
"## Dependencies: \n",
"*Disclaimer: This tutorial has been tested on Windows and MacOS using the specifications identified below.* \n",
"+ #### Python Version 3.8 \n",
" + `h5py` \n",
" + `pyproj` \n",
" + `matplotlib` \n",
" + `pandas` \n",
" + `pyresample` \n",
" + `pandas` \n",
" + `scipy` \n",
" + `gdal` \n",
" + `jupyter notebook`\n",
"\n",
"- A [NASA Earthdata Login](https://urs.earthdata.nasa.gov/) account is required to download the data used in this tutorial. You can create an account at the link provided. \n",
"- The [Geospatial Data Abstraction Library](http://www.gdal.org/) (GDAL) is required.\n",
"\n",
"---\n",
"# Procedures: \n",
"## Getting Started: \n",
"#### 1. This tutorial uses data from ECOSTRESS Version 1, including an ECO3ETPTJPL.001 (and accompanying ECO1BGEO.001) observation from August 05, 2018. You can download the files directly from the [LP DAAC Data Pool](https://e4ftl01.cr.usgs.gov/ECOSTRESS/) at:\n",
" - [ECOSTRESS_L1B_GEO_00468_007_20180805T220314_0601_03.h5](https://e4ftl01.cr.usgs.gov/ECOSTRESS/ECO1BGEO.001/2018.08.05/ECOSTRESS_L1B_GEO_00468_007_20180805T220314_0601_03.h5) \n",
" - [ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04.h5](https://e4ftl01.cr.usgs.gov/ECOSTRESS/ECO3ETPTJPL.001/2018.08.05/ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04.h5) \n",
"**A [NASA Earthdata Login](https://urs.earthdata.nasa.gov/) account is required to download the data used in this tutorial. You can create an account at the link provided.** \n",
"\n",
" - Ancillary Files Needed: \n",
" - The AmeriFlux Latent Heat data used in Section 4 can be downloaded via a [csv file](https://git.earthdata.nasa.gov/projects/LPDUR/repos/tutorial-ecostress/raw/tower_data.csv?at=refs%2Fheads%2Fmain). \n",
"\n",
"The `tower_data.csv` file will need to be downloaded into the same directory as the tutorial in order to execute the tutorial. \n",
"#### 2.\tCopy/clone/download the [ECOSTRESS Tutorial repo](https://git.earthdata.nasa.gov/rest/api/latest/projects/LPDUR/repos/tutorial-ecostress/archive?format=zip), or the desired tutorial from the LP DAAC Data User Resources Repository: \n",
" - [Working with ECOSTRESS Evapotranspiration Data in Python Jupyter Notebook](https://git.earthdata.nasa.gov/projects/LPDUR/repos/tutorial-ecostress/browse/ECOSTRESS_Tutorial.ipynb) \n",
"
\n",
"NOTE: This tutorial was developed specifically for the ECOSTRESS Evapotranspiration PT-JPL Level 3, Version 1 HDF5 files and will need to be adapted to work with other ECOSTRESS products.
\n",
"\n",
"***\n",
"## 1a. Import Packages \n",
"#### Import the python packages required to complete this tutorial."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Import packages\n",
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"import h5py\n",
"import os\n",
"from os.path import join\n",
"import pyproj\n",
"import numpy as np\n",
"import pandas as pd\n",
"import datetime\n",
"from dateutil import parser\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.colors import LinearSegmentedColormap\n",
"from pyresample import geometry as geom\n",
"from pyresample import kd_tree as kdt\n",
"from osgeo import gdal, gdal_array, gdalconst, osr\n",
"\n",
"# Set plots inside of the Jupyter Notebook\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 1b. Set Up the Working Environment\n",
"#### The input directory is defined as the current working directory. Note that you will need to have the jupyter notebook and example data (.h5 and .csv) stored in this directory in order to execute the tutorial successfully."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"input directory:\n",
"D:/Downloads/ecoworkshop/\n",
"output directory:\n",
"D:\\Downloads\\ecoworkshop\\output\\\n"
]
}
],
"source": [
"# Current working directory will be set as the input directory\n",
"inDir = \"D:/Downloads/ecoworkshop/\" #os.getcwd() + os.sep \n",
"print(\"input directory:\\n{}\".format(inDir))\n",
"\n",
"# Set output directory\n",
"outDir = os.path.normpath(os.path.split(inDir)[0] + os.sep + 'output') + os.sep \n",
"print(\"output directory:\\n{}\".format(outDir))\n",
"\n",
"# Create output directory\n",
"if not os.path.exists(outDir): \n",
" os.makedirs(outDir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Make sure that the ECOSTRESS .h5 data files, and Ameriflux ET data file (.csv) are located in the input directory printed above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 1c. Retrieve Files\n",
"#### Make sure that the ECO1BGEO and ECO3ETPTJPL .h5 files listed in the directions have been downloaded to the `inDir` defined above to follow along in the tutorial."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"os.chdir(inDir)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"geolocation:\n",
"ECOSTRESS_L1B_GEO_00468_007_20180805T220314_0601_03.h5\n",
"products:\n",
"ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04.h5\n"
]
}
],
"source": [
"# List directory contents and create lists of ECOSTRESS HDF5 files (GEO, ET)\n",
"geoList = [file for file in os.listdir() if file.endswith('.h5') and 'GEO' in file]\n",
"print(\"geolocation:\\n{}\".format(\"\\n\".join(geoList)))\n",
"ecoList = [file for file in os.listdir() if file.endswith('.h5') and 'GEO' not in file]\n",
"print(\"products:\\n{}\".format(\"\\n\".join(ecoList)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### The standard format for ECOSTRESS filenames is as follows:\n",
"> **ECOSTRESS_L3_ET_PT-JPL**: Product Type \n",
"**00468**: Orbit number; starting at start of mission, ascending equatorial crossing \n",
"**007**: Scene ID; starting at first scene of first orbit \n",
"**20180805T220314**: Date and time of data start: YYYYMMDDThhmmss \n",
"**0601**: Build ID of software that generated product, Major+Minor (2+2 digits) \n",
"**04**: Product version number (2 digits) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"
2. Importing and Interpreting Data
\n",
"\n",
"***\n",
"## 2a. Open an ECOSTRESS HDF5 File\n",
"#### Read in an ECOSTRESS HDF5 file using the `h5py` package."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04\n"
]
}
],
"source": [
"f = h5py.File(ecoList[0]) # Read in ECOSTRESS HDF5 file\n",
"ecoName = ecoList[0].split('.h5')[0] # Keep original filename\n",
"print(ecoName)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 2b. Subset SDS Layers and read SDS Metadata\n",
"#### Identify and generate a list of all the SDS layers in the HDF5 file."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EVAPOTRANSPIRATION PT-JPL/ETcanopy\n",
"EVAPOTRANSPIRATION PT-JPL/ETdaily\n",
"EVAPOTRANSPIRATION PT-JPL/ETinst\n",
"EVAPOTRANSPIRATION PT-JPL/ETinstUncertainty\n",
"EVAPOTRANSPIRATION PT-JPL/ETinterception\n",
"EVAPOTRANSPIRATION PT-JPL/ETsoil\n",
"L3_ET_PT-JPL Metadata/AncillaryFileAerosolOpticalDepth\n",
"L3_ET_PT-JPL Metadata/AncillaryFileAirTemperatureNWP\n",
"L3_ET_PT-JPL Metadata/AncillaryFileAirTemperatureRS\n",
"L3_ET_PT-JPL Metadata/AncillaryFileAlbedo\n"
]
}
],
"source": [
"# Create a list of all SDS inside of the .h5 file\n",
"eco_objs = []\n",
"f.visit(eco_objs.append)\n",
"ecoSDS = [str(obj) for obj in eco_objs if isinstance(f[obj], h5py.Dataset)] \n",
"for dataset in ecoSDS[0:10]: \n",
" print(dataset)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Below, subset the SDS list to the two layers needed for comparison with the ground-based AmeriFlux data, `ETinst` and `ETinstUncertainty`."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ETinst\n",
"ETinstUncertainty\n"
]
}
],
"source": [
"# Subset list to ETinst and ETinstUncertainty\n",
"sds = ['ETinst', 'ETinstUncertainty']\n",
"ecoSDS = [dataset for dataset in ecoSDS if dataset.endswith(tuple(sds))]\n",
"for dataset in ecoSDS:\n",
" print(dataset.split('/')[-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"
3. Performing Swath2grid Conversion
\n",
"\n",
"#### Resample the native ECOSTRESS swath data to a grid with defined coordinate reference system (CRS). \n",
"***\n",
"## 3a. Import Geolocation File\n",
"#### The `latitude` and `longitude` arrays from the ECO1BGEO product for the same ECOSTRESS orbit/scene ID are needed to perform the swath2grid conversion on the ECO3ETPT-JPL file. "
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ECOSTRESS_L1B_GEO_00468_007_20180805T220314_0601_03.h5\n",
"ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04.h5\n"
]
}
],
"source": [
"# Find the matching ECO1BGEO file from the file list\n",
"geo = [geoFile for geoFile in geoList if ecoList[0].split('ECOSTRESS_L3_ET_PT-JPL_')[-1].split('T')[0] in geoFile]\n",
"print(geo[0])\n",
"print(ecoList[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Read in the ECO1BGEO file, search for the `latitude` and `longitude` SDS, and import into Python as arrays."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(5632, 5400)\n"
]
}
],
"source": [
"# Open Geo File\n",
"g = h5py.File(geo[0])\n",
"geo_objs = []\n",
"g.visit(geo_objs.append)\n",
"\n",
"# Search for lat/lon SDS inside data file\n",
"latSD = [str(obj) for obj in geo_objs if isinstance(g[obj], h5py.Dataset) and '/latitude' in obj]\n",
"lonSD = [str(obj) for obj in geo_objs if isinstance(g[obj], h5py.Dataset) and '/longitude' in obj]\n",
"\n",
"# Open SDS as arrays\n",
"lat = g[latSD[0]][()].astype(float)\n",
"lon = g[lonSD[0]][()].astype(float)\n",
"\n",
"# Read the array dimensions\n",
"dims = lat.shape\n",
"print(dims)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 3b. Define Projection and Output Grid\n",
"#### The `latitude` and `longitude` arrays from the ECO1BGEO product for the same ECOSTRESS orbit/scene ID are needed to perform the swath2grid conversion on the ECO3ETPT-JPL file. \n",
"#### The following sections use the [pyresample package](https://pyresample.readthedocs.io/en/stable/) to resample the ECOSTRESS swath dataset to a grid using nearest neighbor method. This process begins by defining the swath dimensions using the lat/lon arrays below."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[(-118.7239829569036, 34.221086430090274),\n",
" (-121.86703190925778, 36.84456587689783),\n",
" (-118.81973784761588, 39.05327678830593),\n",
" (-115.71697211227038, 36.35603480980171)]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Set swath definition from lat/lon arrays\n",
"swathDef = geom.SwathDefinition(lons=lon, lats=lat)\n",
"swathDef.corners"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Define the coordinates in the middle of the swath, which are used to calculate an estimate of the output rows/columns for the gridded output. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(36.63127602811011, -118.88736212064671)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Define the lat/lon for the middle of the swath\n",
"mid = [int(lat.shape[1] / 2) - 1, int(lat.shape[0] / 2) - 1]\n",
"midLat, midLon = lat[mid[0]][mid[1]], lon[mid[0]][mid[1]]\n",
"midLat, midLon"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Below, `pyproj.Proj` is used to perform a cartographic transformation by defining an Azimuthal Equidistant projection centered on the midpoint of the swath. Once the projection is defined, convert the lower left and upper right corners of the lat/lon arrays to a location (in meters) in the new projection. Lastly, measure the distance between the corners and divide by 70 (meters), the nominal pixel size that we are aiming for. Azimuthal Equidistant projection was chosen here based on the following characteristics of this projection:\n",
"- Units in meters (necessary for defining 70 m pixels) \n",
"- Distances between all points are proportionally correct from center point\n",
"- Azimuth (direction) are correct from the center point"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Define AEQD projection centered at swath center\n",
"epsgConvert = pyproj.Proj(\"+proj=aeqd +lat_0={} +lon_0={}\".format(midLat, midLon))\n",
"\n",
"# Use info from AEQD projection bbox to calculate output cols/rows/pixel size\n",
"llLon, llLat = epsgConvert(np.min(lon), np.min(lat), inverse=False)\n",
"urLon, urLat = epsgConvert(np.max(lon), np.max(lat), inverse=False)\n",
"areaExtent = (llLon, llLat, urLon, urLat)\n",
"cols = int(round((areaExtent[2] - areaExtent[0]) / 70)) # 70 m pixel size\n",
"rows = int(round((areaExtent[3] - areaExtent[1]) / 70))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Use number of rows and columns generated above from the AEQD projection to set a representative number of rows and columns in the Geographic area definition, which will then be translated to degrees below, then take the smaller of the two pixel dims to determine output size and ensure square pixels."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Define Geographic projection\n",
"epsg, proj, pName = '4326', 'longlat', 'Geographic'\n",
"\n",
"# Define bounding box of swath\n",
"llLon, llLat, urLon, urLat = np.min(lon), np.min(lat), np.max(lon), np.max(lat)\n",
"areaExtent = (llLon, llLat, urLon, urLat)\n",
"\n",
"# Create area definition with estimated number of columns and rows\n",
"projDict = pyproj.CRS(\"epsg:4326\")\n",
"areaDef = geom.AreaDefinition(epsg, pName, proj, projDict, cols, rows, areaExtent)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Below, square the pixels by setting the pixel size to the smaller of the x any y values output by the `AreaDefinition`, then use the pixel size to recalculate the number of output cols/rows."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Square pixels and calculate output cols/rows\n",
"ps = np.min([areaDef.pixel_size_x, areaDef.pixel_size_y])\n",
"cols = int(round((areaExtent[2] - areaExtent[0]) / ps))\n",
"rows = int(round((areaExtent[3] - areaExtent[1]) / ps))\n",
"\n",
"# Set up a new Geographic area definition with the refined cols/rows\n",
"areaDef = geom.AreaDefinition(epsg, pName, proj, projDict, cols, rows, areaExtent)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Below, use `pyresample kd_tree`'s [get_neighbour_info](https://pyresample.readthedocs.io/en/latest/swath.html#resampling-from-neighbour-info) to create arrays with information on the nearest neighbor to each grid point. \n",
"#### This is the most computationally heavy task in the swath2grid conversion and using `get_neighbour_info` speeds up the process if you plan to resample multiple SDS within an ECOSTRESS product (compute once instead of for every SDS). "
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# Get arrays with information about the nearest neighbor to each grid point \n",
"index, outdex, indexArr, distArr = kdt.get_neighbour_info(swathDef, areaDef, 210, neighbours=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Above, the function is comparing the swath and area definitions to locate the nearest neighbor (neighbours=1). 210 is the `radius_of_influence`, or the radius used to search for the nearest neighboring pixel in the swath (in meters). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 3c. Read SDS Metadata\n",
"#### List the attributes for the `ETinst` layer, which can then be used to define the fill value and scale factor. "
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"_FillValue = nan\n",
"add_offset = 0.0\n",
"long_name = Instantaneous Latent Heat Flux\n",
"scale_factor = 1.0\n",
"units = W/m^2\n",
"valid_max = 2000.0\n",
"valid_min = 0.0\n"
]
}
],
"source": [
"# Read in ETinst and print out SDS attributes\n",
"s = ecoSDS[0]\n",
"ecoSD = f[s][()]\n",
"for attr in f[s].attrs:\n",
" if type(f[s].attrs[attr]) == np.ndarray:\n",
" print(f'{attr} = {f[s].attrs[attr][0]}')\n",
" else:\n",
" print(f'{attr} = {f[s].attrs[attr].decode(\"utf-8\")}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Extract the scale factor, add offset and fill value from the SDS metadata."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"nan"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f[s].attrs['_FillValue'][0]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"# Read SDS attributes and define fill value, add offset, and scale factor if available\n",
"try:\n",
" fv = int(f[s].attrs['_FillValue'])\n",
"except KeyError:\n",
" fv = None\n",
"except ValueError:\n",
" fv = f[s].attrs['_FillValue'][0]\n",
"try:\n",
" sf = f[s].attrs['_Scale'][0]\n",
"except:\n",
" sf = 1\n",
"try:\n",
" add_off = f[s].attrs['_Offset'][0]\n",
"except:\n",
" add_off = 0\n",
"try:\n",
" units = f[s].attrs['units'].decode(\"utf-8\")\n",
"except:\n",
" units = 'none'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 3d. Perform K-D Tree Resampling\n",
"#### Remember that the resampling has been split into two steps. In section 3b. arrays containing the nearest neighbor to each grid point were created. The second step is to use those arrays to retrieve a resampled result. "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# Perform K-D Tree nearest neighbor resampling (swath 2 grid conversion)\n",
"ETgeo = kdt.get_sample_from_neighbour_info('nn', areaDef.shape, ecoSD, index, outdex, indexArr, fill_value=None)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Above, resample the swath `ecoSD` array using nearest neighbor (already calculated in section 3b. and defined above as the `index`, `outdex`, and `indexArr`), and also set the fill value that was defined in section 3c. \n",
"#### Below, define the geotransform for the output (upper left x, horizontal pixel size, rotation, upper left y, rotation, vertical pixel size)."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[-121.86703190925778,\n",
" 0.0006302582963630697,\n",
" 0,\n",
" 39.05327678830593,\n",
" 0,\n",
" -0.0006302582963630697]"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Define the geotransform \n",
"gt = [areaDef.area_extent[0], ps, 0, areaDef.area_extent[3], 0, -ps]\n",
"gt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 3e. Basic Image Processing\n",
"#### Apply the scale factor and add offset and set the fill value defined in the previous section on the resampled data."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"ETgeo = ETgeo * sf + add_off # Apply Scale Factor and Add Offset\n",
"ETgeo[ETgeo == fv * sf + add_off] = fv # Set Fill Value"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Rerun steps 3c - 3e for `ETinstUncertainty`."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"s = ecoSDS[1]\n",
"ecoSD = f[s][()]\n",
"try:\n",
" fv = int(f[s].attrs['_FillValue'])\n",
"except KeyError:\n",
" fv = None\n",
"except ValueError:\n",
" fv = f[s].attrs['_FillValue'][0]\n",
"try:\n",
" sf = f[s].attrs['_Scale'][0]\n",
"except:\n",
" sf = 1\n",
"try:\n",
" add_off = f[s].attrs['_Offset'][0]\n",
"except:\n",
" add_off = 0\n",
"UNgeo = kdt.get_sample_from_neighbour_info('nn', areaDef.shape, ecoSD, index, outdex, indexArr, fill_value=None)\n",
"UNgeo = UNgeo * sf + add_off\n",
"UNgeo[UNgeo == fv * sf + add_off] = fv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"
4. Exporting Results
\n",
"\n",
"***\n",
"## 4a. Set Up a Dictionary\n",
"#### In this section, create a dictionary containing each of the arrays that will be exported as GeoTIFFs."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"# Set up dictionary of arrays to export\n",
"outFiles = {'ETinst': ETgeo, 'ETinstUncertainty': UNgeo}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"## 4b. Define CRS and Export as GeoTIFFs\n",
"#### Now that the data have been imported and resampled into a gridded raster array, export the results as GeoTIFFs using a `for` loop in this section."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"fv = np.nan"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"output file:\n",
"D:\\Downloads\\ecoworkshop\\output\\ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04_ETinst.tif\n",
"\n",
"output file:\n",
"D:\\Downloads\\ecoworkshop\\output\\ECOSTRESS_L3_ET_PT-JPL_00468_007_20180805T220314_0601_04_ETinstUncertainty.tif\n",
"\n"
]
}
],
"source": [
"# Loop through each item in dictionary created above\n",
"for file in outFiles:\n",
" \n",
" # Set up output name\n",
" outName = join(outDir, '{}_{}.tif'.format(ecoName, file))\n",
" print(\"output file:\\n{}\\n\".format(outName))\n",
" \n",
" # Get driver, specify dimensions, define and set output geotransform\n",
" height, width = outFiles[file].shape\n",
" driv = gdal.GetDriverByName('GTiff')\n",
" dataType = gdal_array.NumericTypeCodeToGDALTypeCode(outFiles[file].dtype)\n",
" d = driv.Create(outName, width, height, 1, dataType)\n",
" d.SetGeoTransform(gt)\n",
" \n",
" # Create and set output projection, write output array data\n",
" # Define target SRS\n",
" srs = osr.SpatialReference()\n",
" srs.ImportFromEPSG(int(epsg))\n",
" d.SetProjection(srs.ExportToWkt())\n",
" srs.ExportToWkt()\n",
" \n",
" # Write array to band\n",
" band = d.GetRasterBand(1)\n",
" band.WriteArray(outFiles[file])\n",
" \n",
" # Define fill value if it exists, if not, set to mask fill value\n",
" if fv is not None and fv != 'NaN':\n",
" band.SetNoDataValue(fv)\n",
" else:\n",
" try:\n",
" band.SetNoDataValue(outFiles[file].fill_value)\n",
" except:\n",
" pass\n",
" band.FlushCache()\n",
" d, band = None, None "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"
5. Combining ECOSTRESS and AmeriFlux Tower Data
\n",
"\n",
"***\n",
"## 5a. Loading Tables with Pandas\n",
"#### In this section, begin by highlighting how to open a csv file using the `pandas` package. \n",
"> The AmeriFlux tower data was provided by Mike Goulden for the AmeriFlux US-CZ3 tower. The csv includes half-hourly observations of Latent Heat (W/m$^{2}$) for the same day as the ECOSTRESS observation. "
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"# Import csv file with AmeriFlux data and drop NaNs\n",
"towerData = pd.read_csv('tower_data.csv')\n",
"towerData = towerData.dropna()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Next, use the `parser` package and a lambda function to go through each time stamp and reformat to date and time objects."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"