<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Analysis Ready Data Tutorial Part 1
Step1:  Step 1
Step2:  Step 2
Step3:  As we can see, the footprints (rectangles) do not exactly match the AOI. Indeed, none of them cover the AOI. We don't care about pixels outside of the AOI, so we are going to want to clip the imagery to the AOI (to remove pixels outside the AOI).
Step 3
Step4:  Step 3.1
Step5:  Step 3.2
Step6:  Option 2
Step7:  Step 4
Step8:  Step 4.2
Step9:  Step 4.3
Step10:  Step 5
Step13:  5.2
<ASSISTANT_TASK:>
Python Code:
import datetime
import json
import os
from pathlib import Path
from pprint import pprint
import time
from zipfile import ZipFile
import numpy as np
from planet import api
from planet.api import filters
import rasterio
from rasterio import plot
from shapely.geometry import MultiPolygon, shape
Explanation: Analysis Ready Data Tutorial Part 1: Introduction and Best Practices
Time-series analysis (e.g. change detection and trend detection) is a powerful application of satellite imagery. However, a great deal of processing is required to prepare imagery for analysis. Analysis Ready Data (ARD), preprocessed time-series stacks of overhead imagery, allow for time-series analysis without any additional processing of the imagery. See Analysis Data Defined for an excellent introduction and discussion on ARD.
This tutorial shows how Planet APIs can simplify production of ARD by demonstrating best practices and then by walking through a real world use case. This tutorial is targeted to users who have little to no geospatial knowledge but have experience working with APIs. The goal of this tutorial is to teach the user the how and whys of using the Data and Orders APIs to create and interpret ARD for both use cases. This first part of the totorial focuses on best practices. The following part will focus on a real-world use case.
APIs
The two Planet APIs that are used in creation of ARD are the Data API and the Orders API. The Data API is used to search for imagery according to certain criteria based on the use case, and the Orders API is used to process that imagery into ARD that can be fed directly into time-series analysis.
Data API
The first step in using the Data API is identifying the search criteria. This is specifying answers to the following questions regarding the use case time-series analysis:
* What is the time range?
* What product item type is desired?
* What is the area of interest (geographic region)?
* What percentage of pixels need to be usable?
* etc.
While time range is likely pretty trivial to determine, product item, area of interest, and usable pixels may take a little bit of work. Let's dive into each further
Product Item Type
The product item type refers to the sensor source (aka satellite type) and basic processing desired. This decision is highly dependent on application, as coverage, revisit rate, spectral bands, and resolution differ between products. A good overview of the products available in the Planet Data API is provided on the Planet Imagery and Archive page (look for the link to product specs for details). For most frequent revisit rate, we will use the PS satellite. Experience has shown that customers most often use the scene (vs the orthotile) product. Therefore, this tutorial will focus on the PSScene product.
Area of Interest
The area of interest is the geographic region for the analysis, given as GeoJSON. If you are familiar with JSON, the format of GeoJSON will likely be easy to grasp. According to <geojson.org>, "GeoJSON is a format for encoding a variety of geographic data structures." The specific geographic data structure we are interested is a Polygon. The GeoJSON wikipedia page gives some great examples of the data structures for the various GeoJSON data structures. 
Some care needs to be given to describing the position for each coordinate in a GeoJSON geometry. Each position is specified as (longitude, latitude) or (easting, northing) (order really matters here!). Also, this may be a surprise, but the same point on earth can have different (longitude, latitude) values based on the spatial reference system. Basically, a spatial reference system describes where something is in the real world. But there are thousands of different spatial reference systems. These are separated into two categories: geographic and projected. Geographic coordinate systems model the earth as an ellipsoid (this model is called the datum) and describe positions on the surface of the earth (coordinates) in terms of the prime meridian and angle of measure. Projected coordinate systems take this a step further by projecting the geographic coordinates (quite three-dimensional) into two dimensions (the projection). Different projections preserve different properties, such as area, angles, or direction for north. There is a rich area of discovery, discussion, and even a little teasing (thanks xkcd!) in the world of spatial reference systems.
GeoJSON only supports one spatial reference system, WGS84. This is a geographic coordinate system describing locations in latitude, longitude. However, many web mapping applications use the Web Mercator projected coordinate system to to describe locations. Confusingly, Web Mercator also describes locations in latitude, longitude. But a (longitude, latitude) GeoJSON position given in Web Mercator will not end up where you expect if it is not first projected into WGS84. The easiest way to define an aoi that is described in WGS84 is to draw it in Planet Explorer and click the little button "Download AOI as GeoJSON."
Usable Pixels
Taking pictures of the earth from space ain't easy. There are a lot of steps and a lot of things have to go right to get a clear shot. Therefore, not every pixel in every image taken from space is useful for analysis. For one, images taken from space are projected into a spatial reference system (discussed briefly above), a process that introduces some NoData (aka outside of the image footprint) pixels into the resulting image. Additionally, clouds cover a great deal of the earth and create cloudy pixels when imaged. While some applications can use cloudy pixels, others cannot. Therefore, the type of pixels that are determined to be 'usable' are often application-specific. To support definition of usable pixels, and filtering based on that definition, Planet provides Usable Data Masks along with Usable Data entries in the imagery metadata. For more details on the Usable Data Mask, check out Clear for Analysis with Planet’s New Usable Data Masks. You can also find great examples for working with Usable Data Masks in the UDM2 notebooks. For more information on the Usable Data metadata entries, see Usable Data in Planet imagery.
Orders API
The core decision around using the orders api is which product bundle to use. This is the starting point for all processing and there are a lot of options. Once the product is determined, the processing steps (aka tools and toolchains) are defined. Finally, the logistics of the delivery of the imagery are ironed out.
Product Bundle
To enable time-series analysis, ARD imagery must be processed so that imagery is consistant across days, months, and possibly years. This means correcting for differences in camera sensitivities, the relative location of the sun, and the atmospheric conditions. The Analytic Radiance (analytic) product bundle provides imagery corrected for difference in camera sensitivities and location of the sun (as radiance) and can also remove the effect of the sun's spectrum by applying the reflectanceCoefficient value given in the imagery metadata. However, the Analytic Surface Reflectance (analytic_sr) product bundle removes the effect of atmospheric conditions while also converting to reflectance. Therefore, the Analytic Surface Reflectance is the ideal product for ARD and the one we will use here.
Tools and Toolchains
The Tools and Toolchains functionality in the Orders API are the key to seamlessly creating ARD. Through the API, one can define the pre-processing steps for the data before it is delivered. Given proper definition of the tools and toolchains, data that is delivered is analysis-ready.
Delivery
There are a few options for delivery that cater to different use cases. Imagery can be downloaded directly or delivered to cloud storage. When imagery is downloaded, the user can poll for when the order is ready or notifications (e-mail or webhooks) can be used. Additionally, the imagery can be delivered as a zip archive.
Best Practices
Now that we have a basic understanding of the Data and Orders APIs, let's put them to use creating ARD. This will be a demonstration of best practices.
For this tutorial we will use the planet python client (documentation). This client simplifies interactions with the various Planet APIs and also includes a command-line interface.
The first step in using the planet python client is initializing the client with the user API key. Each user (or organization) has their own unique API key. Information on finding the API key and running Docker so that the API key is available in the notebooks is given in this repository's README.
The next step is building functionality for searching the Data API with the client. This consists of building the search query and then running the search. In building this functionality, we will use test information for data ranges and AOIs to test the functionality, but we will build the functions so those pieces can be changed by the end user. Fo this step we will demonstrate use of  both the python api and the cli.
The third step is building functionality for processing and delivery with the Orders API. At the time of the creation of this tutorial (June 2019), the Orders API functionality best supported (e.g. documented) in the command-line interface (CLI). Therefore, we will use the client CLI for this portion of the tutorial.
Finally, once we have downloaded the order, we will unzip it and visualize the resulting imagery.
To summarize, these are the steps:
1. Initialize API client
1. Search Data API
1. Submit Order
1. Download Order
1. Unzip and View Order
Import Dependencies
End of explanation
# if your Planet API Key is not set as an environment variable, you can paste it below
API_KEY = os.environ.get('PL_API_KEY', 'PASTE_YOUR_KEY_HERE')
client = api.ClientV1(api_key=API_KEY)
Explanation: Step 1: Initialize API client
End of explanation
# define test data for the filter
test_start_date = datetime.datetime(year=2019,month=4,day=1)
test_stop_date = datetime.datetime(year=2019,month=5,day=1)
# iowa crops aoi
test_aoi_geom = {
    "type": "Polygon",
    "coordinates": [
        [
            [-93.299129, 42.699599],
            [-93.299674, 42.812757],
            [-93.288436, 42.861921],
            [-93.265332, 42.924817],
            [-92.993873, 42.925124],
            [-92.993888, 42.773637],
            [-92.998396, 42.754529],
            [-93.019154, 42.699988],
            [-93.299129, 42.699599]
        ]
    ]
}
# create an api request from the search specifications
def build_request(aoi_geom, start_date, stop_date):
    '''build a data api search request for clear PSScene 4-Band imagery'''
    item_type = 'PSScene'
    query = filters.and_filter(
        filters.geom_filter(aoi_geom),
        filters.range_filter('clear_percent', gte=90),
        filters.date_range('acquired', gt=start_date),
        filters.date_range('acquired', lt=stop_date)
    )
    return filters.build_search_request(query, ['PSScene'])
request = build_request(test_aoi_geom, test_start_date, test_stop_date)
print(request)
# search the data api
def search_data_api(request, client, limit=500):
    result = client.quick_search(request)
    
    # this returns a generator
    return result.items_iter(limit=limit)
items = list(search_data_api(request, client))
print(len(items))
# check out an item just for fun
# pprint(items[0])
# visualize a scene footprint
footprints = [shape(i['geometry']) for i in items]
footprints[0]
# visualize subset of footprints and aoi
MultiPolygon([shape(test_aoi_geom), *footprints[:5]])
Explanation: Step 2: Search Data API
The goal of this step is to get the scene ids that meet the search criteria for this use case.
End of explanation
# work with just a subset of the items in the interest of bandwidth
test_items = items[:2]
# filter to item ids
ids = [i['id'] for i in test_items]
ids
# specify the PSScene 4-Band surface reflectance product
# make sure to get the *_udm2 bundle so you get the udm2 product
# note: capitalization really matters in item_type when using planet client orders api
item_type = 'PSScene'
bundle = 'analytic_sr_udm2'
Explanation: As we can see, the footprints (rectangles) do not exactly match the AOI. Indeed, none of them cover the AOI. We don't care about pixels outside of the AOI, so we are going to want to clip the imagery to the AOI (to remove pixels outside the AOI).
Step 3: Submit Order
Now that we have the scene ids, we can create the order. The output of this step is a single zip file that contains all of the scenes that meet our criteria.
Because this is a demo, we don't download all of the scene ids.
End of explanation
# specify tools
# clip to AOI
clip_tool = {'clip': {'aoi': test_aoi_geom}}
# convert to NDVI
bandmath_tool = {'bandmath': {
    "pixel_type": "32R",
    "b1": "(b4 - b3) / (b4+b3)"
}}
tools = [clip_tool, bandmath_tool]
pprint(tools)
# specify a name
name = 'tutorial_order'
Explanation: Step 3.1: Build Orders Toolchain
End of explanation
orders_request = {
    'name': name,
    'products': [{
        'item_ids': ids,
        'item_type': item_type,
        'product_bundle': bundle
    }],
    'tools': tools,
    'delivery': {
        'single_archive': True,
        'archive_filename':'{{name}}_{{order_id}}.zip',
        'archive_type':'zip'
    },
        'notifications': {
                   'email': False
    },
}
# pprint(orders_request, indent=4)
order_info = client.create_order(orders_request).get()
order_id = order_info['id']
order_id
Explanation: Step 3.2: Submit Order
Option 1: Client Python API
The first option for submitting an order is using the planet client python api.
End of explanation
# zip up entire order into one file
ziptype = 'order'
# format the ids for use with the CLI
cli_ids = ','.join([i['id'] for i in test_items])
# save tools definition to file
tools_file = 'tools.json'
with open(tools_file, 'w') as dst:
    dst.write(json.dumps(tools))
    
order_info_file = 'order.json'
# submit the order and save the response to a file so we can get the order id
!set -x;planet orders create \
    --id $cli_ids \
    --item-type $item_type \
    --bundle $bundle \
    --zip $ziptype \
    --tools $tools_file \
    --name $name | tee $order_info_file
# read the order id
with open(order_info_file, 'r') as src:
    order_info = json.load(src)
    
order_id = order_info['id']
order_id
Explanation: Option 2: Client CLI
In some instances, using the CLI to submit orders may be desired.
End of explanation
def poll_for_success(order_id, client, num_loops=50):
    count = 0
    while(count < num_loops):
        count += 1
        order_info = client.get_individual_order(order_id).get()
        state = order_info['state']
        print(state)
        success_states = ['success', 'partial']
        if state == 'failed':
            raise Exception(response)
        elif state in success_states:
            break
        
        time.sleep(10)
        
poll_for_success(order_id, client)
Explanation: Step 4: Download Order
To download the order from the orders api, we will use the planet python client CLI. It would be nice to use the python client python api for this step but, as of the writing of this tutorial, support for the orders api in the planet client python api has been confusing. The CLI download status output is also confusing and possibly wrong but the CLI does easily and successfully download the order.
When we download an order, we always get a manifest.json file. Therefore, while we only ordered one file (an order zip), we will download two files. The manifest is very useful and we will use it to locate the zip file we ordered and downloaded.
Step 4.1: Wait Until Order is Successful
Before we can download the order, the order has to be prepared on the server. Though it would be ideal if the CLI download functionality took care of this, right now it will error out if the order hasn't already been run successfully. So, for now, we will do our own polling and waiting until the order is successful.
End of explanation
demo_data_dir = os.path.join('data', 'demo')
# make the download directory if it doesn't exist
Path(demo_data_dir).mkdir(parents=True, exist_ok=True)
!set -x;planet orders download --dest $demo_data_dir $order_id
Explanation: Step 4.2: Run Download
For this step we will use the orders CLI because it is the easiest and best supported way to download via the planet client for now.
One thing to watch out for: if you have already downloaded an order, it won't be available for re-download. You will have to resubmit the order and get a new order id. What is tricky is that the manifest is still downloaded, but the other files are not.
End of explanation
!ls data/demo
def get_download_locations(download_dir):
    manifest_file = os.path.join(download_dir, 'manifest.json')
    
    with open(manifest_file, 'r') as src:
        manifest = json.load(src)
    
    # uncomment to see the manifest
#     pprint(manifest)
        
    locations = [os.path.join(download_dir, f['path'])
                 for f in manifest['files']]
    return locations
locations = get_download_locations(demo_data_dir)
pprint(locations)
# lets just double check to see if the zip file got downloaded
!ls data/demo
Explanation: Step 4.3: Get Downloaded File Location(s)
We use the downloaded order manifest to find the downloaded file locations. The manifest is saved in the download directory.
End of explanation
def unzip(filename):
    location = Path(filename)
    
    zipdir = location.parent / location.stem
    with ZipFile(location) as myzip:
        myzip.extractall(zipdir)
    return zipdir
zipdir = unzip(locations[0])
zipdir
def get_unzipped_files(zipdir):
    filedir = zipdir / 'files'
    filenames = os.listdir(filedir)
    return [filedir / f for f in filenames]
file_paths = get_unzipped_files(zipdir)
pprint(file_paths)
Explanation: Step 5: Unzip and View Order
In this step we will simply unzip the order and view the downloaded images and their usable data masks.
5.1: Unzip Order
We will unzip the order into a directory named after the file, then we will find the downloaded files (they are in a files subdirectory)
End of explanation
def get_image_and_udm_files(file_paths):
    files = [str(p) for p in file_paths]
    
    # the image files are tiffs and are identified with '_SR_' in the name
    img_id = '_AnalyticMS_SR_'
    imgfiles = [f for f in files
                if f.endswith('.tif') and img_id in f]
    
    # get associated udm files for image files
    # each image has a unique id at the beginning of the name
    imgroots = [str(f).split(img_id)[0] for f in imgfiles]
    
    # the udm files are identified with '_udm2' in the name
    udmfiles = [next(f for f in files if f.startswith(r + '_udm2'))
                for r in imgroots]
    
    return imgfiles, udmfiles
imgfiles, udmfiles = get_image_and_udm_files(file_paths)
pprint(imgfiles)
pprint(udmfiles)
# read UDM2 file
def read_notclear(udm2_filename):
    with rasterio.open(udm2_filename) as img:
        # the first band is the clear/not clear band
        mask=img.read(1)
        not_clear = mask == 0
        return not_clear
    
udmfile = udmfiles[0]
not_clear = read_notclear(udmfile)
# there is an issue where some udms aren't the same size as the images
# to deal with this just cut off any trailing rows/columns
# this isn't ideal as it can result in up to one pixel shift in x or y direction
def crop(img, shape):
    return img[:shape[0], :shape[1]]
def read_ndvi(img_filename, not_clear):
    with rasterio.open(imgfile) as img:
        # ndvi is a single-band image
        band = img.read(1)
        
        # crop image and mask to same size
        img_shape = min(band.shape, not_clear.shape)
        ndvi = np.ma.array(crop(band, img_shape), mask=crop(not_clear, img_shape))
    return ndvi
    
imgfile = imgfiles[0]
ndvi = read_ndvi(imgfile, not_clear)
# set up NDVI visualization
# copied from:  https://stackoverflow.com/a/48598564
import matplotlib.pyplot as plt
import matplotlib.colors as colors
The NDVI values will range from -1 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/
    Credit: https://stackoverflow.com/a/48598564
    
    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):    
        # Note that I'm ignoring clipping and other edge cases here.
        result, is_scalar = self.process_value(value)
        x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
        return np.ma.array(np.interp(value, x, y), mask=result.mask, copy=False)
def show_ndvi(ndvi):
    fig = plt.figure(figsize=(20,10))
    ax = fig.add_subplot(111)
    # diverging color scheme chosen from https://matplotlib.org/users/colormaps.html
    cmap = plt.cm.RdYlGn 
    mmin = np.nanmin(ndvi)
    mmax = np.nanmax(ndvi)
#     print((mmin, mmax))
    mid = 0
    cax = ax.imshow(ndvi, cmap=cmap, clim=(mmin, mmax),
                    norm=MidpointNormalize(midpoint=mid,vmin=mmin, vmax=mmax))
    ax.axis('off')
    ax.set_title('Normalized Difference Vegetation Index', fontsize=18, fontweight='bold')
    cbar = fig.colorbar(cax, orientation='horizontal', shrink=0.65)
    plt.show()
    
# show_ndvi(ndvi)
for imgfile, udmfile in zip(imgfiles, udmfiles):
    show_ndvi(read_ndvi(imgfile, read_notclear(udmfile)))
Explanation: 5.2: Visualize Images
In this section we will find the image files and their associated UDMs and we will visualize them.
The first band of the UDM2 file is the clear/not-clear band. 0: not-clear, 1: clear.
End of explanation
<END_TASK>