<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Let's pull it all together to do something cool.
Let's reuse a lot of our code to make a movie of our travel around Portland.
We'll first select a bunch of recent scenes, activate, and download them.
After that we'll create a mosaic, a path, and trace the path through the moasic. 
We'll use the path to crop subregions, save them as images, and create a video.
First step is to trace our AOI and a path through it.
Step1:  Query the API
Now we'll save the geometry for our AOI and the path.
We'll also filter and cleanup our data just like before.
Step2:  Just like before we clean up our data and distill it down to just the scenes we want.
Step3:  To make sure we are good we'll visually inspect the scenes in our slippy map.
Step12:  This is from the previous notebook. We are just activating and downloading scenes.
Step13:  Perform the actual activation ... go get coffee
Step14:  Downloand the scenes
Step15:  Now, just like before, we will mosaic those scenes.
It is easier to call out using subprocess and use the command line util.
Just iterate through the files and drop them into a single file portland_mosaic.tif
Step16:  Let's take a look at what we got
Step17:  Now we are going to write a quick crop function.
this function takes in a, scene, a center position, and the width and height of a window.
We'll use numpy slice notation to make the crop.
Let's pick a spot and see what we get.
Step19:  Now to figure out how our lat/long values map to pixels.
The next thing we need is a way to map from a lat and long in our slippy map to the pixel position in our image. 
We'll use what we know about the lat/long of the corners of our image to do that. 
We'll ask GDAL to tell us the extents of our scene and the geotransofrm.
We'll then apply the GeoTransform from GDAL to the coordinates that are the extents of our scene. 
Now we have the corners of our scene in Lat/Long
Step20:  Here we'll call the functions we wrote.
First we open the scene and get the width and height.
Then from the geotransorm we'll reproject those points to lat and long.
Step21:  Now we'll do a bit of hack.
That bit above is precise but complext, we are going to make everything easier to think about. 
We are going to linearize our scene, which isn't perfect, but good enough for our application.
What this function does is take in a given lat,long, the size of the image, and the extents as lat,lon coordinates.
For a given pixel we map it's x and y values to the value between a given lat and long and return the results.
Now we can ask, for a given lat,long pair what is the corresponding pixel.
Step22:  Let's check our work
First we'll create a draw point function that just puts a red dot at given pixel.
We'll get our scene, and map all of the lat/long points in our path to pixel values.
Finally we'll load our image, plot the points and show our results
Step23:  Now things get interesting....
Our path is just a few waypoint but to make a video we need just about every point between our waypoints.
To get all of the points between our waypoints we'll have to write a little interpolation script. 
Interpolation is just a fancy word for nicely space points bewteen or waypoints, we'll call the space between each point as our "velocity."
If we were really slick we could define a heading vector and and build a spline so the camera faces the direction of heading. Our approach is fine as the top of the frame is always North, which makes reckoning easy.
Once we have our interpolation function all we need to do is to crop our large mosaic at each point in our interpolation point list and save it in a sequential file.
Step24:  Before we generate our video frames, let's check our work
We'll load our image. 
Build the interpolated waypoints list.
Draw the points on the image using our draw_point method.
Plot the results
Step25:  Now let's re-load the image and run the scene maker.
Step26:  Finally, let's make a movie.
Our friend AVConv, which is like ffmpeg is a handy command line util for transcoding video.
AVConv can also convert a series of images into a video and vice versa.
We'll set up our command and use subprocess to make the call.
<ASSISTANT_TASK:>
Python Code:
# Basemap Mosaic (v1 API)
mosaicsSeries = 'global_quarterly_2017q1_mosaic'
# Planet tile server base URL (Planet Explorer Mosaics Tiles)
mosaicsTilesURL_base = 'https://tiles0.planet.com/experimental/mosaics/planet-tiles/' + mosaicsSeries + '/gmap/{z}/{x}/{y}.png'
# Planet tile server url
mosaicsTilesURL = mosaicsTilesURL_base + '?api_key=' + api_keys["PLANET_API_KEY"]
# Map Settings 
# Define colors
colors = {'blue': "#009da5"}
# Define initial map center lat/long
center = [45.5231, -122.6765]
# Define initial map zoom level
zoom = 11
# Set Map Tiles URL
planetMapTiles = TileLayer(url= mosaicsTilesURL)
# Create the map
m = Map(
    center=center, 
    zoom=zoom,
    default_tiles = planetMapTiles # Uncomment to use Planet.com basemap
)
# Define the draw tool type options
polygon = {'shapeOptions': {'color': colors['blue']}}
rectangle = {'shapeOptions': {'color': colors['blue']}} 
# Create the draw controls
# @see https://github.com/ellisonbg/ipyleaflet/blob/master/ipyleaflet/leaflet.py#L293
dc = DrawControl(
    polygon = polygon,
    rectangle = rectangle
)
# Initialize an action counter variable
actionCount = 0
AOIs = {}
# Register the draw controls handler
def handle_draw(self, action, geo_json):
    # Increment the action counter
    global actionCount
    actionCount += 1
    # Remove the `style` property from the GeoJSON
    geo_json['properties'] = {}
    # Convert geo_json output to a string and prettify (indent & replace ' with ")
    geojsonStr = json.dumps(geo_json, indent=2).replace("'", '"')
    AOIs[actionCount] = json.loads(geojsonStr)
    
# Attach the draw handler to the draw controls `on_draw` event
dc.on_draw(handle_draw)
m.add_control(dc)
m
Explanation: Let's pull it all together to do something cool.
Let's reuse a lot of our code to make a movie of our travel around Portland.
We'll first select a bunch of recent scenes, activate, and download them.
After that we'll create a mosaic, a path, and trace the path through the moasic. 
We'll use the path to crop subregions, save them as images, and create a video.
First step is to trace our AOI and a path through it.
End of explanation
areaAOI = AOIs[1]["geometry"]
pathAOI = AOIs[2]["geometry"]
aoi_file ="portland.geojson" 
with open(aoi_file,"w") as f:
    f.write(json.dumps(areaAOI))
# build a query using the AOI and
# a cloud_cover filter that excludes 'cloud free' scenes
old = datetime.datetime(year=2017,month=1,day=1)
query = filters.and_filter(
    filters.geom_filter(areaAOI),
    filters.range_filter('cloud_cover', lt=5),
    filters.date_range('acquired', gt=old)
)
# build a request for only PlanetScope imagery
request = filters.build_search_request(
    query, item_types=['PSScene3Band']
)
# if you don't have an API key configured, this will raise an exception
result = client.quick_search(request)
scenes = []
planet_map = {}
for item in result.items_iter(limit=500):
    planet_map[item['id']]=item
    props = item['properties']
    props["id"] = item['id']
    props["geometry"] = item["geometry"]
    props["thumbnail"] = item["_links"]["thumbnail"]
    scenes.append(props)
scenes = pd.DataFrame(data=scenes)
display(scenes)
print len(scenes)
Explanation: Query the API
Now we'll save the geometry for our AOI and the path.
We'll also filter and cleanup our data just like before.
End of explanation
# now let's clean up the datetime stuff
# make a shapely shape from our aoi
portland = shape(areaAOI)
footprints = []
overlaps = []
# go through the geometry from our api call, convert to a shape and calculate overlap area.
# also save the shape for safe keeping
for footprint in scenes["geometry"].tolist():
    s = shape(footprint)
    footprints.append(s)
    overlap = 100.0*(portland.intersection(s).area / portland.area)
    overlaps.append(overlap)
# take our lists and add them back to our dataframe
scenes['overlap'] = pd.Series(overlaps, index=scenes.index)
scenes['footprint'] = pd.Series(footprints, index=scenes.index)
# now make sure pandas knows about our date/time columns.
scenes["acquired"] = pd.to_datetime(scenes["acquired"])
scenes["published"] = pd.to_datetime(scenes["published"])
scenes["updated"] = pd.to_datetime(scenes["updated"])
scenes.head()
# Now let's get it down to just good, recent, clear scenes
clear = scenes['cloud_cover']<0.4
good = scenes['quality_category']=="standard"
recent = scenes["acquired"] > datetime.date(year=2017,month=1,day=1)
partial_coverage = scenes["overlap"] > 10
good_scenes = scenes[(good&clear&recent&partial_coverage)]
print good_scenes
Explanation: Just like before we clean up our data and distill it down to just the scenes we want.
End of explanation
# first create a list of colors
colors = ["#ff0000","#00ff00","#0000ff","#ffff00","#ff00ff","#00ffff","#ff0000","#00ff00","#0000ff","#ffff00","#ff00ff","#00ffff"]
# grab our scenes from the geometry/footprint geojson
# Chane this number as needed
footprints = good_scenes[0:10]["geometry"].tolist()
# for each footprint/color combo
for footprint,color in zip(footprints,colors):
    # create the leaflet object
    feat = {'geometry':footprint,"properties":{
            'style':{'color': color,'fillColor': color,'fillOpacity': 0.2,'weight': 1}},
            'type':u"Feature"}
    # convert to geojson
    gjson = GeoJSON(data=feat)
    # add it our map
    m.add_layer(gjson)
# now we will draw our original AOI on top 
feat = {'geometry':areaAOI,"properties":{
            'style':{'color': "#FFFFFF",'fillColor': "#FFFFFF",'fillOpacity': 0.5,'weight': 1}},
            'type':u"Feature"}
gjson = GeoJSON(data=feat)
m.add_layer(gjson)   
m 
Explanation: To make sure we are good we'll visually inspect the scenes in our slippy map.
End of explanation
def get_products(client, scene_id, asset_type='PSScene3Band'):    
    
    Ask the client to return the available products for a 
    given scene and asset type. Returns a list of product 
    strings
    
    out = client.get_assets_by_id(asset_type,scene_id)
    temp = out.get()
    return temp.keys()
def activate_product(client, scene_id, asset_type="PSScene3Band",product="analytic"):
    
    Activate a product given a scene, an asset type, and a product.
    
    On success return the return value of the API call and an activation object
    
    temp = client.get_assets_by_id(asset_type,scene_id)  
    products = temp.get()
    if( product in products.keys() ):
        return client.activate(products[product]),products[product]
    else:
        return None 
def download_and_save(client,product):
    
    Given a client and a product activation object download the asset. 
    This will save the tiff file in the local directory and return its 
    file name. 
    
    out = client.download(product)
    fp = out.get_body()
    fp.write()
    return fp.name
def scenes_are_active(scene_list):
    
    Check if all of the resources in a given list of
    scene activation objects is read for downloading.
    
    return True
    retVal = True
    for scene in scene_list:
        if scene["status"] != "active":
            print "{} is not ready.".format(scene)
            return False
    return True
def load_image4(filename):
    Return a 4D (r, g, b, nir) numpy array with the data in the specified TIFF filename.
    path = os.path.abspath(os.path.join('./', filename))
    if os.path.exists(path):
        with rasterio.open(path) as src:
            b, g, r, nir = src.read()
            return np.dstack([r, g, b, nir])
        
def load_image3(filename):
    Return a 3D (r, g, b) numpy array with the data in the specified TIFF filename.
    path = os.path.abspath(os.path.join('./', filename))
    if os.path.exists(path):
        with rasterio.open(path) as src:
            b,g,r,mask = src.read()
            return np.dstack([b, g, r])
        
def get_mask(filename):
    Return a 1D mask numpy array with the data in the specified TIFF filename.
    path = os.path.abspath(os.path.join('./', filename))
    if os.path.exists(path):
        with rasterio.open(path) as src:
            b,g,r,mask = src.read()
            return np.dstack([mask])
def rgbir_to_rgb(img_4band):
    Convert an RGBIR image to RGB
    return img_4band[:,:,:3]
Explanation: This is from the previous notebook. We are just activating and downloading scenes.
End of explanation
to_get = good_scenes["id"][0:10].tolist()
to_get = sorted(to_get)
activated = []
# for each scene to get
for scene in to_get:
    # get the product 
    product_types = get_products(client,scene)
    for p in product_types:
        # if there is a visual productfor p in labels:
        if p == "visual": # p == "basic_analytic_dn"
            print "Activating {0} for scene {1}".format(p,scene)
            # activate the product
            _,product = activate_product(client,scene,product=p)
            activated.append(product)
Explanation: Perform the actual activation ... go get coffee
End of explanation
tiff_files = []
asset_type = "_3B_Visual"
# check if our scenes have been activated
if scenes_are_active(activated):
    for to_download,name in zip(activated,to_get):
        # create the product name
        name = name + asset_type + ".tif"
        # if the product exists locally
        if( os.path.isfile(name) ):
            # do nothing 
            print "We have scene {0} already, skipping...".format(name)
            tiff_files.append(name)
        elif to_download["status"] == "active":
            # otherwise download the product
            print "Downloading {0}....".format(name)
            fname = download_and_save(client,to_download)
            tiff_files.append(fname)
            print "Download done."
        else:
            print "Could not download, still activating"
else:
    print "Scenes aren't ready yet"
print tiff_files 
Explanation: Downloand the scenes
End of explanation
subprocess.call(["rm","portland_mosaic.tif"])
commands = ["gdalwarp", # t
           "-t_srs","EPSG:3857",
           "-cutline",aoi_file,
           "-crop_to_cutline",
           "-tap",
            "-tr", "3", "3"
           "-overwrite"]
output_mosaic = "portland_mosaic.tif"
for tiff in tiff_files:
    commands.append(tiff)
commands.append(output_mosaic)
print " ".join(commands)
subprocess.call(commands)
Explanation: Now, just like before, we will mosaic those scenes.
It is easier to call out using subprocess and use the command line util.
Just iterate through the files and drop them into a single file portland_mosaic.tif
End of explanation
merged = load_image3(output_mosaic)
plt.figure(0,figsize=(18,18))
plt.imshow(merged)
plt.title("merged")
Explanation: Let's take a look at what we got
End of explanation
def crop_to_area(scene,x_c,y_c,w,h):
    tlx = x_c-(w/2)
    tly = y_c-(h/2)
    brx = x_c+(w/2)
    bry = y_c+(h/2)
    return scene[tly:bry,tlx:brx,:]
plt.figure(0,figsize=(3,4))
plt.imshow(crop_to_area(merged,3000,3000,640,480))
plt.title("merged")
Explanation: Now we are going to write a quick crop function.
this function takes in a, scene, a center position, and the width and height of a window.
We'll use numpy slice notation to make the crop.
Let's pick a spot and see what we get.
End of explanation
# Liberally borrowed from this example
# https://gis.stackexchange.com/questions/57834/how-to-get-raster-corner-coordinates-using-python-gdal-bindings
def GetExtent(gt,cols,rows):
    
    Get the list of corners in our output image in the format 
    [[x,y],[x,y],[x,y]]
    
    ext=[]
    # for the corners of the image
    xarr=[0,cols]
    yarr=[0,rows]
    for px in xarr:
        for py in yarr:
            # apply the geo coordiante transform 
            # using the affine transform we got from GDAL
            x=gt[0]+(px*gt[1])+(py*gt[2])
            y=gt[3]+(px*gt[4])+(py*gt[5])
            ext.append([x,y])
        yarr.reverse()
    return ext
def ReprojectCoords(coords,src_srs,tgt_srs):
    trans_coords=[]
    # create a transform object from the source and target ref system
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    for x,y in coords:
        # transform the points
        x,y,z = transform.TransformPoint(x,y)
        # add it to the list. 
        trans_coords.append([x,y])
    return trans_coords
Explanation: Now to figure out how our lat/long values map to pixels.
The next thing we need is a way to map from a lat and long in our slippy map to the pixel position in our image. 
We'll use what we know about the lat/long of the corners of our image to do that. 
We'll ask GDAL to tell us the extents of our scene and the geotransofrm.
We'll then apply the GeoTransform from GDAL to the coordinates that are the extents of our scene. 
Now we have the corners of our scene in Lat/Long
End of explanation
# TLDR: pixels => UTM coordiantes => Lat Long 
raster=output_mosaic
# Load the GDAL File
ds=gdal.Open(raster)
# get the geotransform
gt=ds.GetGeoTransform()
# get the width and height of our image
cols = ds.RasterXSize
rows = ds.RasterYSize
# Generate the coordinates of our image in utm
ext=GetExtent(gt,cols,rows)
# get the spatial referencec object 
src_srs=osr.SpatialReference()
# get the data that will allow us to move from UTM to Lat Lon. 
src_srs.ImportFromWkt(ds.GetProjection())
tgt_srs = src_srs.CloneGeogCS()
extents = ReprojectCoords(ext,src_srs,tgt_srs)
print extents
Explanation: Here we'll call the functions we wrote.
First we open the scene and get the width and height.
Then from the geotransorm we'll reproject those points to lat and long.
End of explanation
def poor_mans_lat_lon_2_pix(lon,lat,w,h,extents):
    # split up our lat and longs 
    lats = [e[1] for e in extents]
    lons = [e[0] for e in extents]
    # calculate our scene extents  max and min
    lat_max = np.max(lats)
    lat_min = np.min(lats)    
    lon_max = np.max(lons)
    lon_min = np.min(lons) 
    # calculate the difference between our start point
    # and our minimum
    lat_diff = lat-lat_min
    lon_diff = lon-lon_min
    # create the linearization
    lat_r = float(h)/(lat_max-lat_min)
    lon_r = float(w)/(lon_max-lon_min) 
    # generate the results. 
    return int(lat_r*lat_diff),int(lon_r*lon_diff)
Explanation: Now we'll do a bit of hack.
That bit above is precise but complext, we are going to make everything easier to think about. 
We are going to linearize our scene, which isn't perfect, but good enough for our application.
What this function does is take in a given lat,long, the size of the image, and the extents as lat,lon coordinates.
For a given pixel we map it's x and y values to the value between a given lat and long and return the results.
Now we can ask, for a given lat,long pair what is the corresponding pixel.
End of explanation
def draw_point(x,y,img,t=40):
    h,w,d = img.shape
    y = h-y
    img[(y-t):(y+t),(x-t):(x+t),:] = [255,0,0]
h,w,c = merged.shape
waypoints = [poor_mans_lat_lon_2_pix(point[0],point[1],w,h,extents) for point in pathAOI["coordinates"]]
print waypoints
merged = load_image3(output_mosaic)
[draw_point(pt[1],pt[0],merged) for pt in waypoints]
plt.figure(0,figsize=(18,18))
plt.imshow(merged)
plt.title("merged")
Explanation: Let's check our work
First we'll create a draw point function that just puts a red dot at given pixel.
We'll get our scene, and map all of the lat/long points in our path to pixel values.
Finally we'll load our image, plot the points and show our results
End of explanation
def interpolate_waypoints(waypoints,velocity=10.0):
    retVal = []
    last_pt = waypoints[0]
    # for each point in our waypoints except the first
    for next_pt in waypoints[1:]:
        # calculate distance between the points
        distance = np.sqrt((last_pt[0]-next_pt[0])**2+(last_pt[1]-next_pt[1])**2)
        # use our velocity to calculate the number steps.
        steps = np.ceil(distance/velocity)
        # linearly space points between the two points on our line
        xs = np.array(np.linspace(last_pt[0],next_pt[0],steps),dtype='int64')
        ys = np.array(np.linspace(last_pt[1],next_pt[1],steps),dtype='int64')
        # zip the points together
        retVal += zip(xs,ys)
        # move to the next point
        last_pt = next_pt
    return retVal
def build_scenes(src,waypoints,window=[640,480],path="./movie/"):
    count = 0 
    # Use opencv to change the color space of our image.
    src = cv2.cvtColor(src, cv2.COLOR_BGR2RGB)
    # define half our sampling window. 
    w2 = window[0]/2
    h2 = window[1]/2
    # for our source image get the width and height
    h,w,d = src.shape
    for pt in waypoints:
        # for each point crop the area out.
        # the y value of our scene is upside down. 
        temp = crop_to_area(src,pt[1],h-pt[0],window[0],window[1])
        # If we happen to hit the border of the scene, just skip
        if temp.shape[0]*temp.shape[1]== 0:
            # if we have an issue, just keep plugging along
            continue
        # Resample the image a bit, this just makes things look nice. 
        temp = cv2.resize(temp, (int(window[0]*0.75), int(window[1]*.75))) 
        # create a file name
        fname = os.path.abspath(path+"img{num:06d}.png".format(num=count))
        # Save it
        cv2.imwrite(fname,temp)
        count += 1
Explanation: Now things get interesting....
Our path is just a few waypoint but to make a video we need just about every point between our waypoints.
To get all of the points between our waypoints we'll have to write a little interpolation script. 
Interpolation is just a fancy word for nicely space points bewteen or waypoints, we'll call the space between each point as our "velocity."
If we were really slick we could define a heading vector and and build a spline so the camera faces the direction of heading. Our approach is fine as the top of the frame is always North, which makes reckoning easy.
Once we have our interpolation function all we need to do is to crop our large mosaic at each point in our interpolation point list and save it in a sequential file.
End of explanation
# load the image
merged = load_image3(output_mosaic)
# interpolate the waypoints
interp = interpolate_waypoints(waypoints)
# draw them on our scene
[draw_point(pt[1],pt[0],merged) for pt in interp]
# display the scene
plt.figure(0,figsize=(18,18))
plt.imshow(merged)
plt.title("merged")
Explanation: Before we generate our video frames, let's check our work
We'll load our image. 
Build the interpolated waypoints list.
Draw the points on the image using our draw_point method.
Plot the results
End of explanation
os.system("rm ./movie/*.png")
merged = load_image3(output_mosaic)
build_scenes(merged,interp)
Explanation: Now let's re-load the image and run the scene maker.
End of explanation
# avconv -framerate 30 -f image2 -i ./movie/img%06d.png -b 65536k out.mpg
os.system("rm ./movie/*.png")
framerate = 30
output = "out.mpg"
command = ["avconv","-framerate", str(framerate), "-f", "image2", "-i", "./movie/img%06d.png", "-b", "65536k", output]
os.system(" ".join(command))
Explanation: Finally, let's make a movie.
Our friend AVConv, which is like ffmpeg is a handy command line util for transcoding video.
AVConv can also convert a series of images into a video and vice versa.
We'll set up our command and use subprocess to make the call.
End of explanation
<END_TASK>