Planet's analytic satellite data products are reported in units of radiance: that is to say, every pixel in an analytic-type PlanetScope image answers the question, "How much light was captured over this spot of ground?" The mathematical representation of this radiance is $W*m^{-2}*sr^{-1}$.
But over the course of a day or year, the number of photons that the sun shines on the scene rises and falls. If you naively compare the radiance values of two scenes over the same spot on Earth from the same satellite but a few weeks (or even just hours) apart, you will likely find dramatic radiance differences even if nothing on the ground has changed.
In addition to this variation, each of Planet Labs' 150+ satllites have small amounts of variation in their spectral filters which yields slight differences in radiance measurements, even from two satellites taking pictures of the same exact place at the same exact moment.
To correct all this radiance variation you would have to do a lot of math using the exact location and local time of day to find the angle to the sun and sun-earth distance, then compute a solar irradiance model to estimate how many photons of each band are present in the image, and finally convolve that spectrum with the spectral response of the individual satellite to yield the number of photons of each band that are actually recorded by that sensor. Dividing by this number normalizes the measured brightness to the brightness of the sun at that time and place through that particular filter, yielding a much more comparable number: reflectance.
One type of reflectance correction is known as "Top of Atmosphere Reflectance" (TOAR). TOAR is extremely useful because it is an apples-to-apples comparable number from any satellite over any location that does not change with time of day or time of year unless the content of the scene changes. It is very commonly used in geospatial applications that compute spectral indices such as NDVI. If you'd like to learn how to compute the NDVI index, see this in-class exercise.
TOAR is so broadly useful that Planet has already done the math for us, and provides all the coefficients necessary to convert a radiance image into a reflectance image.
The four coefficients $a, b, c, d$ are calculated and provided with every analytic image that Planet provides and can be used as simple scaling factors to convert from radiance to reflectance. Their values change with the image's local time of day and time of year, and do so uniquely per satellite.
In this exercise, you'll learn perform a basic Radiance to Reflectance calculation on PlanetScope imagery using Python. Here are the steps to follow:
$PL_API_KEY.item-type: PSScene4Band; asset-type: analyticFor this exercise you'll need a 4-band PlanetScope Analytioc product. You can search for & download your own data, or use the demo data provided in-class. If you choose to use the demo data, skip to Step 2.
To search for your own data, you'll first need to define an Area of Interest (AOI). http://geojson.io is a free browser-based tool that makes generating a GeoJSON-formatted AOI easy.
Once that's done, use one of the following methods to search for & download data:
With all of the above, you'll want to filter for 4-Band PlanetScope data (item type: PSScene4Band) and download the associated analytic product (asset type: analytic).
In [ ]:
# To use Planet's CLI from this Notebook, begin your line as follows:
!planet data
# Here is an example of using Planet's CLI to search for a known item id:
# !planet data download --item-type PSScene4Band --asset-type analytic,analytic_xml --string-in id 20170623_180038_0f34
If you prefer to use Planet's API directly via Python, this search & download quickstart Notebook may be useful.
In [ ]:
# To use Planet's API, you'll probably begin by importing your favorite HTTP toolkit, e.g.:
import requests
from requests.auth import HTTPBasicAuth
# Your Planet API key is available in this Notebook as an env variable, e.g.:
import os
PLANET_API_KEY = os.getenv('PL_API_KEY')
If you prefer to browse for images visually, log in to your Planet account and use Planet Explorer to search for PlanetScope imagery. You'll want to make sure to set the Source filter to show only 4-band PlanetScope Scene results.
You can click here for an example search showing 4-band PlanetScope data in California's Central Valley.
Regardless of the path you chose to obtain data for this exercise, once you have successfully aquired a 4-band PlanetScope analytic-type GeoTIFF and associated .xml file, place the files in the data/ directory adjacent to this Notebook. The first of these two files is a GeoTIFF, the image you requested with spatial reference data embedded. The second file is a metadata file for that image that includes the data you'll need to convert from radiance to reflectance.
In this step, you'll use Rasterio, a Python library for reading and writing geospatial raster datasets, to open the raster image you downloaded (the .tif file). Then you'll extract the data from the red and near-infrared bands and load the band data into arrays that you can manipulate using Python's NumPy libary.
Note: in PlanetScope 4-band images, the band order is BGRN: (1) Blue, (2) Green, (3) Red, (4) Near-infrared.
In [ ]:
import rasterio
import numpy as np
filename = "data/20170623_180038_0f34_3B_AnalyticMS.tif"
# Load bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(filename) as src:
band_blue_radiance = src.read(1)
with rasterio.open(filename) as src:
band_green_radiance = src.read(2)
with rasterio.open(filename) as src:
band_red_radiance = src.read(3)
with rasterio.open(filename) as src:
band_nir_radiance = src.read(4)
Before we can convert to reflectance, the conversion coefficients from the metadata file (the .xml file) must be extracted. For the purposes of this exercise, we've included this process here:
In [ ]:
from xml.dom import minidom
xmldoc = minidom.parse("data/20170623_180038_0f34_3B_AnalyticMS_metadata.xml")
nodes = xmldoc.getElementsByTagName("ps:bandSpecificMetadata")
# XML parser refers to bands by numbers 1-4
coeffs = {}
for node in nodes:
bn = node.getElementsByTagName("ps:bandNumber")[0].firstChild.data
if bn in ['1', '2', '3', '4']:
i = int(bn)
value = node.getElementsByTagName("ps:reflectanceCoefficient")[0].firstChild.data
coeffs[i] = float(value)
print("Conversion coefficients:", coeffs)
Note that the coefficients are all of order 1e-5, and that the coefficient for NIR is significantly higher than the coefficient for blue. This is a big deal if your use case involves performing band math because a pixel with a NIR/blue ratio of 1.0 in the radiance image will have a NIR/blue ratio of 3.35/1.929=1.73 in the reflectance image.
Most spectral indices are defined in terms of reflectance, not radiance.
Radiance is measured in SI units: $W/m^2$. Reflectance is a ratio from 0 to 1. The conversion is performed as a per-band scalar multiplication:
In [ ]:
# Multiply the current values in each band by the TOA reflectance coefficients
band_blue_reflectance = band_blue_radiance * coeffs[1]
band_green_reflectance = band_green_radiance * coeffs[2]
band_red_reflectance = band_red_radiance * coeffs[3]
band_nir_reflectance = band_nir_radiance * coeffs[4]
In [ ]:
# check the results by taking a look at the min-max range of the radiance (before) vs reflectance (after)
import numpy as np
# for our sample purposes, we'll check just the red band here
radiance_min, radiance_max = np.amin(band_red_radiance), np.amax(band_red_radiance)
reflectance_min, reflectance_max = np.amin(band_red_reflectance), np.amax(band_red_reflectance)
print("Red band radiance is from {} to {}".format(radiance_min, radiance_max))
print("Red band reflectance is from {} to {}".format(reflectance_min, reflectance_max))
Finally, we'll save the calculated reflectance values to a new image file, making sure the new image file has the same geospatial metadata as the original GeoTIFF we downloaded.
A note:
Reflectance is generally defined as a floating point number between 0 and 1, but image file formats are much more commonly stored as unsigned integers. A common practice in the industry is to multiply the reflectance value by 10,000, then save the result as a file with data type uint16.
In [ ]:
# get the metadata of original GeoTIFF:
meta = src.meta
print(meta)
# set the source metadata as kwargs we'll use to write the new data:
kwargs = meta
# update the 'dtype' value to uint16:
kwargs.update(dtype='uint16')
# Use the min & max values we previously computed for the red band
print("Red band values before scaling: {}-{}".format(reflectance_min, reflectance_max))
# As noted above, scale reflectance value by a factor of 10k:
scale = 10000
blue_ref_scaled = scale * band_blue_reflectance
green_ref_scaled = scale * band_green_reflectance
red_ref_scaled = scale * band_red_reflectance
nir_ref_scaled = scale * band_nir_reflectance
# Compute new min & max values for the scaled red band, just for comparison
red_min_scaled = np.amin(red_ref_scaled)
red_max_scaled = np.amax(red_ref_scaled)
print("Red band values before scaling: {}-{}".format(red_min_scaled, red_max_scaled))
# Convert the type to 'uint16'
from rasterio import uint16
red = red_ref_scaled.astype(uint16)
green = green_ref_scaled.astype(uint16)
blue = blue_ref_scaled.astype(uint16)
nir = nir_ref_scaled.astype(uint16)
# Finally, write band calculations to a new GeoTIFF file
with rasterio.open('reflectance.tif', 'w', **kwargs) as dst:
dst.write_band(1, red)
dst.write_band(2, green)
dst.write_band(3, blue)
dst.write_band(4, nir)
Congratulations! You've learned how to apply coefficient metadata to a GeoTIFF in order to atmospherically correct radiance values, converting them to reflectance.
From here, you could use your new reflectance image to generate indices such as an NDVI. To learn more about how to do that, see band-math-generate-ndvi/generate-ndvi-exercise.ipynb.