Accessing the Planet Labs data API from R

Planet labs

Planet labs provides some amazing data for monitoring earth’s ecosystems. Their constellation of satellites provides panchromatic and multiband images of the entire earth almost daily, with spatial resolution down to 80cm. This is an amazing resource for monitoring ecosystems and global change impacts

They offer a 14-day free trial, and you can apply for free usage if you wish to use the data for academic, or non-profit purposes.

Planet + Shiny

I recently built a Shiny application that used Planet data to examine regions within the fynbos biome that are not behaving as would be expected if the vegetation was healthy. You can read more about this project here, and look at the code here.

Planet data API mechanics

Access to the data is most convinient using their command line or Python tools, but to access the Planet data API from R I had to write some functions mimicing the Python API. In case you need to access the Planet data API from R, or what to include Planet data in a shiny application here is a little description of the function. It provides access to the search endpoint of the Planet data API. This API call will return the itemID of the Planet data that match the search critera. After getting these itemIDs we can either download the binary data for further proceessing using the Planet data API download endpoint, or use the itemIDs to add the Planet data to an interactive map which will allow users to explore the region of interest, or even compare images taken at different dates. I have not written a function to access the data download API endpoint (I will get around to this), but after showing you the function to access itemIDs from the search API endpoint I will show you how to use these to add the Planet imagery to a leaflet map which you can interact with or - as I did - add to a Shiny application. You can read all about the mechanics of the Planet API here.

Here is my function to access the search endpoint:

# bbox = bounding box (shapefile)
# date_end = date of last image (date)
# date_start = date of first image (date)
# cloud_lim = cloud cover limit (double)
# cover_lim = area coverage limit (double)
# item_name = type of image (string)

get_planet <- function(bbox, date_end, date_start, cloud_lim=0.2, cover_lim=0.5, item_name="PSOrthoTile")
  {
    #convert shapefile to geojson
    #shapefile of bounding box must be EPSG:4326 Projection

    geo_json_geometry <- list(
      type=jsonlite::unbox("Polygon"),
      coordinates = list(list(
        c(bbox@xmin,
          bbox@ymin),
        c(bbox@xmin,
          bbox@ymax),
        c(bbox@xmax,
          bbox@ymax),
        c(bbox@xmax,
          bbox@ymin),
        c(bbox@xmin,
          bbox@ymin)
      ))
    )


    # filter for items the overlap with our chosen geometry
    geometry_filter <- list(
      type= jsonlite::unbox("GeometryFilter"),
      field_name= jsonlite::unbox("geometry"),
      config= geo_json_geometry
    )

    #we will search for images for up to a month beforethe date we are interested in

    dategte <- paste0(date_start,"T00:00:00.000Z")
    datelte <- paste0(date_end,"T00:00:00.000Z")

    # filter images by daterange
    date_range_filter <- list(
      type= jsonlite::unbox("DateRangeFilter"),
      field_name= jsonlite::unbox("acquired"),
      config= list(
        gte= jsonlite::unbox(dategte),
        lte= jsonlite::unbox(datelte))
      )

    # filter by cloud cover
    cloud_cover_filter <- list(
      type= jsonlite::unbox("RangeFilter"),
      field_name= jsonlite::unbox("cloud_cover"),
      config = list(
        lte= jsonlite::unbox(cover_lim))
    )

    # filter by coverage of bounding box
    coverage_filter <- list(
      type= jsonlite::unbox("RangeFilter"),
      field_name= unbox("usable_data"),
      config = list(
        gte= jsonlite::unbox(cover_lim))
    )

    # combine filters
    filter_configs <- list(
      type= jsonlite::unbox("AndFilter"),
      config = list(date_range_filter, cloud_cover_filter,geometry_filter, coverage_filter)
    )

    #build request
    search_endpoint_request <- list(
      item_types = item_name,
      filter = filter_configs
    )

    #convert request to JSON
    body_json <- jsonlite::toJSON(search_endpoint_request,pretty=TRUE)

    #API request config
    url <- 'https://api.planet.com/data/v1/quick-search'
    body <- body_json
    api_key <- "myapikey"

    #send API request
    request <- httr::POST(url, body = body, content_type_json(), authenticate(api_key, ""))

    #get request content
    response <- httr::content(request)

    return(response)
}

Planet + leaflet

Let’s see how this works in the wild. I will use Planet data as the basemap for my leaflet map in R. First let’s create some data with which to call the planet API. We need a bounding box and some dates.

library(leaflet)
library(lubridate)
library(httr)
library(jsonlite)
library(XML)
library(dplyr)

#set date limits
date_end <- Sys.Date()
date_start <- today - years(1)

#create a dummy raster and use its extent as bounding box
my_raster<-raster(nrows=100, ncols=100, xmn=20.0, xmx=20.01, ymn=-33.01, ymx=-33.0)
bbox <- extent(my_raster)

#call API
response <- get_planet(bbox, date_end, date_start)

Now that I have a list of items that match my search criteria I can add the data to my interactive map. Often there will be multiple items that match our search criteria. We can access individual items using response$features[[i]] where i is the index of the item we want. Here I will just use the first item and add it to my map.

#get item id of first item
response_id <- response$features[[1]]$id

#create the string to send to XYZ Basemap Tile Service
planet_string <- paste0(
      "https://tiles.planet.com/data/v1/PSOrthoTile/",
      response_id,
      "/{z}/{x}/{y}.png?api_key=myapikey")

#create leaftlet map with planet basemap
my_map = leaflet() %>%
         setView(lng = 20, lat = -33, zoom = 15) %>%
         addTiles(planet_string, attribution = 'Map data @2017 Planet')

my_map