Planet Labs' analytic data products (both Rapideye and Dove) are reported in units of radiance: $W*m^{-2}*sr^{-1}$. That means that every pixel in an analytic tiff has a natural language interpretation: "How much light was captured over this spot of ground?"
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.
Top of Atmosphere Reflectance 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 GIS applications which compute spectral indices such as NDVI or NDWI.
It is so broadly useful that Planet does this math for you 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 guide, you'll perform a basic Radiance to Reflectance calculation on PlanetScope imagery using just a few lines of Python. Here are the steps:
$PL_API_KEY.item-type: PSOrthoTile, REOrthoTile, or PSScene4Band; asset-type: analytic, or basic_analyticFirst, you're going to download a 4-band PlanetScope satellite image of agricultural land in California's Central Valley, captured in mid July 2017 (item-id: 20170623_180038_0f34). You can do this using Planet's Python client to interact with our Data API, or by browsing Planet Explorer, filtering for 4 Band PlanetScope scene (PSScene4Band) or Planetscope ortho tile (PSOrthoTile), and downloading an analytic asset.
Before you download the full image, you can preview a thumbnail of the image via Planet's Data API. (The thumbnails are 256x256 by default, and can be scaled up by passing in a width parameter.)
In [1]:
from IPython.display import Image
Image(url="https://api.planet.com/data/v1/item-types/PSScene4Band/items/20170623_180038_0f34/thumb?width=512")
Out[1]:
Next, you'll use Planet's Python client to download the image. Note: when you run this command, you'll get a stream of messages in your Jupyter notebook as the Python client polls the Data API to determine if the image is activated and ready to download.
In [2]:
!planet data download --item-type PSScene4Band --asset-type analytic,analytic_xml \
--dest data --string-in id 20170623_180038_0f34
Congratulations! You now have two files in your data directory: 20170623_180038_0f34_3B_AnalyticMS.tif and 20170623_180038_0f34_3B_AnalyticMS_metadata.xml. The first file 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 [3]:
import rasterio
import numpy as np
filename = "data/20170623_180038_0f34_3B_AnalyticMS.tif"
# Load red and NIR 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 you can convert to reflectance, you must extract the conversion coefficients from the metadata file you downloaded (the .xml file).
In [4]:
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: {}".format(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 simple per-band scalar multiplication:
In [5]:
# Multiply the Digital Number (DN) 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]
import numpy as np
print("Red band radiance is from {} to {}".format(np.amin(band_red_radiance), np.amax(band_red_radiance)))
print("Red band reflectance is from {} to {}".format(np.amin(band_red_reflectance), np.amax(band_red_reflectance)))
Next, you're going to 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.
In [6]:
# Set spatial characteristics of the output object to mirror the input
kwargs = src.meta
kwargs.update(
dtype=rasterio.uint16,
count = 4)
print("Before Scaling, red band reflectance is from {} to {}".format(np.amin(band_red_reflectance),
np.amax(band_red_reflectance)))
# Here we include a fixed scaling factor. This is common practice.
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
print("After Scaling, red band reflectance is from {} to {}".format(np.amin(red_ref_scaled),
np.amax(red_ref_scaled)))
# Write band calculations to a new raster file
with rasterio.open('data/reflectance.tif', 'w', **kwargs) as dst:
dst.write_band(1, band_blue_reflectance.astype(rasterio.uint16))
dst.write_band(2, band_green_reflectance.astype(rasterio.uint16))
dst.write_band(3, band_red_reflectance.astype(rasterio.uint16))
dst.write_band(4, band_nir_reflectance.astype(rasterio.uint16))
This warrants some explanation: 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 radiance value by 10,000, then save the result as a file with data type uint16.
In the last step, you'll use Matplotlib to visualize the reflectance values you calculated for the PlanetScope scene.
In [8]:
import matplotlib.pyplot as plt
import matplotlib.colors as colors
"""
The reflectance values will range from 0 to 1. You want to use a diverging color scheme to visualize the data,
and you want to center the colorbar at a defined midpoint. The class below allows you to normalize the colorbar.
"""
class MidpointNormalize(colors.Normalize):
"""
Normalise the colorbar so that diverging bars work there way either side from a prescribed midpoint value)
e.g. im=ax1.imshow(array, norm=MidpointNormalize(midpoint=0.,vmin=-100, vmax=100))
Credit: Joe Kington, http://chris35wills.github.io/matplotlib_diverging_colorbar/
"""
def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
self.midpoint = midpoint
colors.Normalize.__init__(self, vmin, vmax, clip)
def __call__(self, value, clip=None):
# I'm ignoring masked values and all kinds of edge cases to make a
# simple example...
x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
return np.ma.masked_array(np.interp(value, x, y), np.isnan(value))
# Set min/max values from reflectance range for image (excluding NAN)
min=np.nanmin(band_nir_reflectance)
max=np.nanmax(band_nir_reflectance)
mid=0.20
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(111)
# diverging color scheme chosen from https://matplotlib.org/users/colormaps.html
# note that appending '_r' to the color scheme name reverses it!
cmap = plt.cm.get_cmap('RdGy_r')
cax = ax.imshow(band_nir_reflectance, cmap=cmap, clim=(min, max), norm=MidpointNormalize(midpoint=mid,vmin=min, vmax=max))
ax.axis('off')
ax.set_title('NIR Reflectance', fontsize=18, fontweight='bold')
cbar = fig.colorbar(cax, orientation='horizontal', shrink=0.65)
fig.savefig("data/ref-fig.png", dpi=200, bbox_inches='tight', pad_inches=0.7)
plt.show()
Looking at the NIR reflectance image, we see that the lakes and ponds appear to reflect almost no infrared light at all, but the healthy fields reflect almost 50% of the incoming infrared light!