Getting Started

Objective: Find a place, search for imagery, and rasterize it.

This tutorial demonstrates a few of the basic APIs in the Descartes Labs Platform. We’ll start by figuring out where we want to look. This will utilize the Places functionality. From there, we’ll briefly cover how to search for products using the Metadata API. Finally, we’ll rasterize the available imagery into a numpy array and display it using their built in methods.

import descarteslabs as dl
import matplotlib.pyplot as plt

So you now have access to a giant archive of imagery. First question is – where do you want to look? You might answer this question in many ways, but one of the ways we can help is by providing mechanisms to find shapes of known places. Our Places API has a find method that does fuzzy-matching searches of places. As an example, let’s try to find Taos, New Mexico (a favorite place for Cartesians to go hiking, biking, camping, and skiing).

# Find potential matches
matches = dl.places.find('new-mexico_taos')
matches
[
  {
    u'bbox': [-106.058364, 36.013014, -105.200117, 36.995841],
    u'id': 102081181,
    u'name': u'Taos',
    u'path': u'continent:north-america_country:united-st...n:new-mexico_district:northwest_county:taos',
    u'placetype': u'county',
    u'slug': u'north-america_united-states_new-mexico_northwest_taos'
  }
]
# The first one looks good to me, so lets make that our area of interest.
aoi = matches[0]
# This area of interest just gives us some basic properties such as bounding boxes.
# To access a GeoJSON Geometry object of that place, we call the `Places.shape` method, in this case
# accessing a low-resolution version of this particular shape.
shape = dl.places.shape(aoi['slug'], geom='low')

If our Places API does not have the boundary you need, or you are working with a unique geometry, you can use your own custom GeoJSON shapes with the rest of our APIs seamlessly.

Searching for available imagery

The Descartes Labs platform pulls together data from over 75 different sources. We call each of these distinct offerings products, and they all have their own unique IDs. Even derived or external data you save into the platform has a product ID. The first parameter when searching for data is a list of product IDs to specify what kind of imagery you want returned. To view products you have access to, you can use the catalog interface which allows you to search our products and displays all their metadata information including details about spectral bands, classes, lifespan, and more. The catalog is an interface for our Metadata API. You can get back all the same information programmatically using the Metadata API.

# Return all the products you have access to
products = dl.metadata.products()
print("I currently have access to {} products. Some of these may not be publically available.".format(len(products)))
I currently have access to 177 products. Some of these may not be publically available.
# Get back all band names based on a product ID
bands = dl.metadata.bands('landsat:LC08:PRE:TOAR')
[band['name'] for band in bands]
[u'coastal-aerosol',
 u'blue',
 u'green',
 u'red',
 u'nir',
 u'swir1',
 u'swir2',
 u'cirrus',
 u'tirs1',
 u'alpha',
 u'bright-mask',
 u'cloud-mask',
 u'qa_cirrus',
 u'qa_cloud',
 u'qa_snow',
 u'qa_water']
# There is metadata associated with each band. In this case, we can tell that the "red" band
# is stored as a UInt16 dataset, has a valid range of [0, 10000] which maps to [0, 1.0] in
# Top-of-atmosphere-reflectance.
band_information['red']
{u'color': u'Red',
 u'dtype': u'UInt16',
 u'name': u'red',
 u'nbits': 14,
 u'nodata': None,
 u'physical_range': [0.0, 1.0],
 u'shortname': u'r',
 u'valid_range': [0, 10000]}

Rasterizing imagery

There are two ways to rasterize an image:

  • dl.raster.raster() : Creates a geo-referenced file and returns it in the response
  • dl.raster.ndarray() : Returns a numpy ndarray object

Both of these methods can take either a single scene key or multiple scene keys. If multiple are supplied, sources will be mosaicked together in order, so the first source will be on the bottom and the last on top. We’ll use the ndarray method in this tutorial.

# Collect the id's for each feature
ids = [f['id'] for f in feature_collection['features']]
# Rasterize the features.
#  * Select red, green, blue, alpha
#  * Scale the incoming data with range [0, 10000] down to [0, 4000] (40% TOAR)
#  * Choose an output type of "Byte" (uint8)
#  * Choose 60m resolution
#  * Apply a cutline of Taos county
arr, meta = dl.raster.ndarray(
    ids,
    bands=['red', 'green', 'blue', 'alpha'],
    scales=[[0,4000], [0, 4000], [0, 4000], None],
    data_type='Byte',
    resolution=60,
    cutline=shape['geometry'],
)

# Note: A value of 1 in the alpha channel signifies where there is valid data.
# We use this throughout the majority of our imagery as a standard way of specifying
# valid or nodata regions. This is particularly helpful if a value of 0 in a particular
# band has meaning, rather than specifying a lack of data.
# We'll use matplotlib to make a quick plot of the image.
import matplotlib.pyplot as plt
%matplotlib inline
plt.figure(figsize=[16,16])
plt.imshow(arr)
<matplotlib.image.AxesImage at 0x7f7ca256d050>
https://cdn.descarteslabs.com/descarteslabs-python/getting_started_files/getting_started_19_1.png
# We can choose other false color band combinations, like
# NIR - SWIR1 - SWIR2
arr, meta = dl.raster.ndarray(
    ids,
    bands=['nir', 'swir1', 'swir2', 'alpha'],
    scales=[[0,4000], [0, 4000], [0, 4000], None],
    cutline=shape['geometry'],
    data_type='Byte',
    resolution=60
)
plt.figure(figsize=[16,16])
plt.imshow(arr)
<matplotlib.image.AxesImage at 0x7f7ca1f21850>
https://cdn.descarteslabs.com/descarteslabs-python/getting_started_files/getting_started_20_1.png
# Or even derived bands like NDVI. Here the alpha channel comes in
# particularly useful as a way to select valid data. Here you may want to use
# some of the band information to scale NDVI properly.

valid_range = band_information['ndvi']['valid_range']
physical_range = band_information['ndvi']['physical_range']
print "%s maps to %s" % (valid_range, physical_range)
arr, meta = dl.raster.ndarray(
    [f['id'] for f in feature_collection['features']],
    bands=['ndvi', 'alpha'],
    scales=[[valid_range[0], valid_range[1], physical_range[0], physical_range[1]], None],
    cutline=shape['geometry'],
    data_type='Float32',
    resolution=60
)
[0, 65535] maps to [-1.0, 1.0]
# Here we can make a numpy masked array using alpha == 0 as a nodata mask.
import numpy as np
mask = arr[:, :, 1] == 0
masked_ndvi = np.ma.masked_array(arr[:, :, 0], mask)
plt.figure(figsize=[16,16])
plt.imshow(masked_ndvi, cmap='BrBG', vmin=0, vmax=0.5)
cb = plt.colorbar()
cb.set_label("NDVI")
https://cdn.descarteslabs.com/descarteslabs-python/getting_started_files/getting_started_22_0.png