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. From there, we’ll briefly cover how to search for products using the Metadata API and what a GeoContext is. Next, we will demonstrate how to search for products and filter imagery using the Scenes API. Finally, we’ll rasterize the available imagery into a numpy array and display it using their built in methods.

In [1]: import descarteslabs as dl

In [2]: import matplotlib.pyplot as plt

In [3]: from descarteslabs.vectors import FeatureCollection, properties as p



So you now have access to a giant archive of imagery. First question is – where do you want to look? In general you’ll have to define a geometry of the area that you’re interested in. There are some ways to help you define that area. One is the website http://geojson.io/. Another one is the vectors service which contains various products with geometries, such as states, cities, counties, and more. Some examples:

• States of the USA
• name: us_states_area
• id: 2c3a5d20b8204023ab07218ae37ee31
• Cities of the USA
• name: us_cities_area
• id: d1349cc2d8854d998aa6da92dc2bd24
• Counties of the USA
• name: us_county_area
• id: 1f360e60a2ae4e3687e7d04eea2faaf

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
In [5]: fc = FeatureCollection('1f360e60a2ae4e3687e7d04eea2faaf')

In [6]: matches = [f for f in fc.filter(properties=(p.county_name == 'Taos')).features()]

In [7]: matches
Out[7]:
[Feature({
'geometry': {
'coordinates': (
(
(-106.058364, 36.296978),
(-106.052349, 36.294994),
(-106.00071, 36.27795),
(-105.999389, 36.277728),
(-105.992675, 36.275444),
(-105.898493, 36.244354),
...
)
,),
'type': 'Polygon'
},
'id': u'35055',
'properties': {
u'aff_geo_code': u'0500000US35055',
u'area_land_meters': 5704119046,
u'area_water_meters': 3383126,
u'county_fips_code': u'055',
u'county_gnis_code': u'00933056',
u'county_name': u'Taos',
u'geo_id': u'35055',
u'legal_area_code': u'06',
u'state_fips_code': u'35'
},
'type': 'Feature'
})]

# Let's make Taos county our area of interest
In [8]: aoi = matches[0]

# And convert it into a GeoJSON Feature object (which contains 'geometry')
In [9]: shape = aoi.geojson


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

In [11]: print("I currently have access to {} products. Some of these may not be publically available.".format(len(products)))
I currently have access to 200 products. Some of these may not be publically available.

# Get back all band names based on a product ID

In [13]: [band['name'] for band in bands]
Out[13]:
[u'coastal-aerosol',
u'blue',
u'green',
u'red',
u'nir',
u'swir1',
u'swir2',
u'cirrus',
u'tirs1',
u'alpha',
u'qa_cirrus',
u'qa_cloud',
u'qa_snow',
u'qa_water']


Lets find some Landsat 8 imagery over our AOI using the Scenes API¶

There are two ways to search for and obtain data: the Scenes API, or a combination of the Metadata API and the Raster API. The Scenes API handles the details of scaling, clipping, resolution, coordinate systems, and no data values through a paradigm we call a GeoContext, while the latter requires more knowledge of the imagery you are requesting. We will focus on the Scenes API for the remainder of this tutorial.

# Create a SceneCollection matching our search criteria and view its summary
In [14]: scenes, ctx = dl.scenes.search(shape['geometry'],
....:                                products="landsat:LC08:01:RT:TOAR",
....:                                start_datetime="2017-11-01",
....:                                end_datetime="2018-07-01",
....:                                limit=500)
....:

In [15]: scenes
Out[15]:
SceneCollection of 60 scenes
* Dates: Nov 08, 2017 to Jun 27, 2018
* Products: landsat:LC08:01:RT:TOAR: 60


Inspect the GeoContext of the SceneCollection¶

When a SceneCollection is created, an associated GeoContext is instantiated which intelligently sets defaults for rasterizing the imagery. This includes: cutline or boundary, spatial resolution, and coordinate reference system.

A GeoContext is immutable, but you can overwrite the default values by creating a new context as demonstrated below.

 # View the context of our SceneCollection
In [16]: ctx
Out[16]:
AOI(geometry=<shapely...97f2df90>,
resolution=15,
crs=u'EPSG:32613',
align_pixels=True,
bounds=(-106.058364, 36.013014, -105.200117, 36.996017),
bounds_crs='EPSG:4326',
shape=None)

# Create a copy of the original context, and set the resolution to 60m.
In [17]: lowres_context = ctx.assign(resolution=60)


The Scenes API uses Shapely under the hood. This means you have accesss to all the Shapely functionality such as plotting the geometry.

# Plot the shape of Taos inline using the geometry method.
In [18]: ctx.geometry
Out[18]: <shapely.geometry.polygon.Polygon at 0x7f1597f2df90>


Display a single image¶

The original Scene Collection has 41 Scenes, though no single Scene has data covering our entire AOI. Before we mosaic the data, let’s take a look at a single scene, accessing its data in the form of an ndarray.

# Pick a Scene that has more than a sliver of data within our AOI.
In [19]: scene = scenes[16]

# A Scene or SceneCollection only contains the Metadata for the imagery
# We have to use another method to get the data back, in this case, the ndarray method
# Note that we are passing in the low resolution GeoContext
arr = scene.ndarray("red green blue", lowres_context)
dl.scenes.display(arr, title="Partial Coverage")


Create a mosaic of the SceneCollection¶

Now that we understand what a single Scene is, let’s take a look at the SceneCollection. We know a single image provides insufficient data for our AOI. The image above is missing data in the upper right hand corner of the shape. The mosaic method will give us complete coverage by overlaying scenes into a single raster. Note that where the images overlap, the pixels from the last image in the list simply “cover up” the others. The method returns a NumPy masked array, where pixels outside our AOI or containing invalid data are masked out.

# Request a mosaic of data containing the red green and blue bands
mosaic = scenes.mosaic("red green blue", lowres_context)
dl.scenes.display(mosaic, title="Taos Mosaic")


Create a cloud free composite¶

The mosaic function simply takes the last valid data value for any given pixel location, sometimes resulting in cloudy or undesirable results as seen above. One way to avoid issues with data in individual scenes is to create a composite: a single image of data aggregated from many scenes. To remove clouds, a simple yet often-effective method is to take the median value of each pixel across many images.

We can do that using SceneCollection.stack and NumPy. The stack method returns a 4D masked ndarray, with each Scene in the SceneCollection as a separate “layer”. Then the compositing is a simple call to np.ma.median.

# Request a stack of data containing the red green and blue bands.
In [20]: stack = scenes.stack("red green blue", lowres_context)

In [21]: print("The stack method returns a {} containing the data for the same {} Scenes returned by the search method above.".format(type(stack),len(stack)))
The stack method returns a <class 'numpy.ma.core.MaskedArray'> containing the data for the same 60 Scenes returned by the search method above.

# We import the NumPy library so we can call the median method on our masked array to produce the final dataset
import numpy as np
composite = np.ma.median(stack, axis=0)
dl.scenes.display(composite, title="Taos Composite")


Of course, you can work with any number of bands available for analysis and visualization. To complete the tutorial, let’s use the same process to display a false color image. Areas that appear red indicate healthy vegetation.

false_color_stack = scenes.stack("nir swir1 swir2", lowres_context)

composite = np.ma.median(false_color_stack, axis=0)
dl.scenes.display(composite, title="Taos Infrared Composite")