Comparing Landsat and PlanetScope Scenes

A common and powerful workflow for analysis is comparing or combining information from multiple sources. This tutorial demonstrates comparing a set of Landsat and PlanetScope scenes.

In this tutorial, we download scenes from the two sensors taken on the same day at the same place, visualize them, resample the Landsat scene to match the pixel resolution and locations of the PlanetScope scene, and perform pixel-by-pixel comparison of their near-infrared (NIR) bands.

Our comparison of the NIR bands is identification of the difference between the scaled NIR bands. We find that most differences between the NIR bands are at high-resolution features and are due to the lower resolution of the Landsat scene. However, we also find a few fields show noticeable difference in the NIR bands. We theorize that these differences may be due to the different spectral responses of the two sensors and the fields containing different crop types than the surrounding fields. These results suggest that an analysis utilizing both Landsat8 and PlanetScope scenes may be able to better differentiate crop types than just using one sensor.

Ideas for extending beyond this notebook:

  1. Mask the scenes using the PlanetScope Unusable Data Mask and Landsat Quality Assurance bands
  2. Convert the PlanetScope scene to reflectance and then do not normalize the bands in the difference calculation

Install Dependencies and Set up Notebook

In [1]:
import os
import shutil
from subprocess import check_output, STDOUT, CalledProcessError
from six.moves import urllib

import numpy as np
import rasterio
import rasterio.enums
import rasterio.vrt
import requests

%matplotlib inline

Download Scenes

In the crossovers notebook, we identified many Landsat/PlanetScope crossovers within the same day that occured between January and August 2017. From that list, we are going to download the following set of scenes:

PSOrthoTile Landsat8L1G

These scenes represent a crossover between Landsat8 and PlanetScope on July 25, 2017.

For the Landsat8L1G scene, we download band 5, which is the NIR band (ref).

To download the scenes, we use the planet CLI because it handles activating, waiting for activation, and downloading the file.

We will save the scenes in the data folder. This folder isn't tracked by git so the downloaded image file will not bloat our git repository.

Download PlanetScope Orthotile Scene

In [2]:
# Create a working directory if it does not already exist.
data_folder = 'data'
if not os.path.isdir(data_folder):

In [3]:
# First test if scene file exists, if not, use the Planet commandline tool to download the image.
# This command assumes a bash shell, available in Unix-based operating systems.
!test -f data/644787_1056721_2017-07-25_0f52_BGRN_Analytic.tif || \
    planet data download \
        --item-type PSOrthoTile \
        --asset-type analytic \
        --string-in id 644787_1056721_2017-07-25_0f52 \
        --dest data

In [4]:
# Store the image filename for processing and make sure the file exists.
pl_filename = os.path.join('data', '644787_1056721_2017-07-25_0f52_BGRN_Analytic.tif')
assert os.path.isfile(pl_filename)


Download Landsat NIR Band

Landsat is distributed with each asset as an individual band. Band 5 is the NIR band. Since we only plan on using the NIR band, we will download just that asset.

NOTE: The command below should work, but currently there is a bug in the planet CLI. Until that is fixed, use the curl command to download the asset from the endpoint directly.

In [5]:
# !planet data download \
#     --item-type Landsat8L1G \
#     --asset-type analytic_b5 \
#     --string-in id LC80430332017206LGN00 \
#     --dest data/

# !ls -l --block-size=M data/

In [6]:
# The long way of downloading an asset. Doesn't use the planet CLI

def get_auth():
    auth = (os.environ['PL_API_KEY'], '')
    return auth

def get_asset_info(asset):
    auth = get_auth()
    item_assets_url = ''
    r = requests.get(item_assets_url, auth=auth)
    assert r.status_code == 200
    assets_resp = r.json()
    return assets_resp[asset]

def activate(asset):
    asset_info = get_asset_info(asset)
    activation_link = asset_info['_links']['activate']
    auth = get_auth()
    r =, auth=auth)
    assert r.status_code == 204

def download(asset, path):
    if not os.path.isfile(path):
        asset_info = get_asset_info(asset)
        download_link = asset_info['location']
        urllib.request.urlretrieve(download_link, path)

asset = 'analytic_b5'

# it may take a while for the asset to be activated. If this fails, wait a while and rerun
landsat_filename = os.path.join('data', 'LC80430332017206LGN00_{}.tif'.format(asset))
download(asset, landsat_filename)

In [7]:
# store the Landsat band filename for processing and make sure the file exists
l8_nir_orig_filename = os.path.join('data', 'LC80430332017206LGN00_analytic_b5.tif')
assert os.path.isfile(l8_nir_orig_filename)


Warp Landsat Scene to Planet Scene

For pixel-by-pixel comparison of the Landsat and Planet scenes, the two scenes must have the same extents, the pixels must be the same size, and the pixels must exactly overlap each other. This can all be achieved by warping one scene to the other scene. We warp the Landsat scene to the Planet scene because the Planet scene is much higher resolution than the Landsat scene and we want to use the highest resolution possible for our comparison.

The GDAL command-line interface, specifically gdalwarp, is the easiest way to warp one image to another image.

In [8]:
def prepare_l8_band(band_filename, dst_filename, out_filename):
    """Project, crop, and resample landsat 8 band to match dst_filename image."""
    # Open the Landsat8 band file.
    with as src:
        # Open the Planet image to use as a template.
        with as pl:
            # A Virtual Raster (VRT) allows for describing raster transformations and delaying computation
            # until data is accessed.  The WarpedVRT() class describes how to warp a dataset, but the
            # reprojection calculations are not transformed, and they are only transformed
            # These parameters define a target image in the same coordinate reference system as the
            # PlanetScope image and the same number of rows and columns.
            warp_params = {
                'transform': pl.transform,
                'width': pl.width,
                'height': pl.height
            with rasterio.vrt.WarpedVRT(src, **warp_params) as vrt:
                # Using the PlanetScope image as a reference, create an output.  The PlanetScope image
                # has 4 bands but the Landsat 8 image only has 1.  Update the output image profile with
                # this information.
                out_profile = pl.profile.copy()
                out_profile['count'] = vrt.count
                with, 'w', **out_profile) as dst:
                    # This '' triggers the warp calculations.
                    data =
                    # Write data to output file.

l8_nir_filename = os.path.join('data', 'LC80430332017206LGN00_analytic_b5_warped.tif')
prepare_l8_band(l8_nir_orig_filename, pl_filename, l8_nir_filename)

Load and Visualize NIR Bands

Before moving on to comparing the NIR bands, we load the NIR bands from each scene and visualize the bands.

Instead of including a huge chunk of code for visualizing the NIR band, we put that code in a local script ( and import here.

Ideally, we would use the PlanetScope Unusable Data Mask (UDM) and Landsat Quality Assurance (QA) bands to create the pixel masks. In the interest of keeping this tutorial focused, we perform a simple and brute-force masking where we mask pixels that are equal to the blackfill NoData value, 0.

In [9]:
import visual

In [10]:
# Display Landsat 8 NIR band.

with, 'r') as src:
    data =
    # Generate a masked array with a mask representing 0 pixels.
    mask = data == 0
    l8_nir_band =, mask=mask)

visual.plot_image(l8_nir_band, title='Landsat8 NIR')

In [11]:
# Display PlanetScope NIR band.

with, 'r') as src:
    data =
    mask = data == 0
    pl_nir_band =, mask=mask)

visual.plot_image(pl_nir_band, title='PSOrthotile NIR')