<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Converting PlanetScope Imagery from Radiance to Reflectance
How to use this exercise
Adjacent to this notebook, you'll find a notebook titled convert-radiance-to-reflectance-key.ipynb
Step1:  Option 2
Step2:  Option 3
Step3:  Step 3. Extract the Coefficients
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
Step4:  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.
Step 4
Step5:  As a quick check of our calculations, let's print the minimum and maximum values in both the original (radiance) values, as well as our new (reflectance) values. 
(HINT
Step6:  Step 5. Save the Reflectance Image
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
<ASSISTANT_TASK:>
Python Code:
# To use Planet's CLI from this Notebook, begin your line as follows:
!planet data
# < add your own code here >
Explanation: Converting PlanetScope Imagery from Radiance to Reflectance
How to use this exercise
Adjacent to this notebook, you'll find a notebook titled convert-radiance-to-reflectance-key.ipynb: that notebook contains the "answers", while the code in this notebook is incomplete. As you work through this notebook, you'll see some codeblocks that contain None values: replace these None instances with your own code until the codeblocks here run successfully. If you get stuck try to use outside resources like documentation, previous work, or even other students here in this class before taking a peek at the Key.
Background
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 $Wm^{-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.
$$
\begin{equation}
\mathbf{img}{reflectance} =  \begin{vmatrix}
a \
b \
c \
d 
\end{vmatrix}
\times \mathbf{img}{radiance}
\end{equation}
$$
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:
Download a PlanetScope image
Extract data from each spectral band
Extract the coefficients
Convert Radiance to Reflectance
Save the Reflectance
Requirements
Python 2.7 or 3+
Planet's Python Client
rasterio
numpy
matplotlib
Planet API Key, stored as environment variable $PL_API_KEY.
Planet 4-Band Imagery with the following specifications: item-type: PSScene; asset-type: ortho_analytic_4b
Step 1. Download a PlanetScope Image
For 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:
- using Planet's Python CLI to interact with Planet's API from the command line
- using Planet's API directly to search and download
- using the Planet Explorer site to visually search for & download data
With all of the above, you'll want to filter for 4-Band PlanetScope data (item type: PSScene) and download the associated analytic product (asset type: ortho_analytic_4b).
Option 1: Searching & Downloading via CLI
If you choose to use Planet's CLI, you might fight these search and download quickstart guides to be useful.
End of explanation
# 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')
# < add your own code here >
Explanation: Option 2: Searching & Downloading via API
If you prefer to use Planet's API directly via Python, this search & download quickstart Notebook may be useful.
End of explanation
import rasterio
import numpy as np
filename = None
# < update above with your own code >
# Load blue, green, red and NIR bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(filename) as src:
    band_blue_radiance = None
# < update above with your own code >
# ...
# < add your own code here >
Explanation: Option 3: Searching & Downloading via Planet Explorer
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.
Success! Data Obtained
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. 
Step 2. Extract the Data from Each Spectral Band
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.
End of explanation
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)
Explanation: Step 3. Extract the Coefficients
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:
End of explanation
# Multiply the current values in each band by the TOA reflectance coefficients
band_blue_reflectance = None
# < update above with your own code >
# ...
# < add your own code here >
Explanation: 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.
Step 4: Convert Radiance to Reflectance
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:
End of explanation
import numpy as np
# < add your own code here >
Explanation: As a quick check of our calculations, let's print the minimum and maximum values in both the original (radiance) values, as well as our new (reflectance) values. 
(HINT: these band values are stored as numpy arrays, so use numpy functions here).
End of explanation
# get the metadata of original GeoTIFF:
# < add your own code here >
# set the source metadata as kwargs we'll use to write the new data:
# < add your own code here >
# update the 'dtype' value to 'uint16':
# < add your own code here >
# As noted above, scale reflectance values by a factor of 10k:
scale = 10000
blue_ref_scaled = None
# < update above with your own code >
# ...
# < add your own code here >
# Compute new min & max values for the scaled red band & print them, just for comparison
# < add your own code here >
# set the type for each band to uint16
# hint: remember, bands are just numpy arrays!
from rasterio import uint16
# < add your own code here >
# Finally, write band calculations to a new raster file 'data/reflectance.tif'
# Hint: rasterio.write_band() to write individual bands
# < add your own code here >
Explanation: Step 5. Save the Reflectance Image
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.
End of explanation
<END_TASK>