Menu Sidebar
Menu

Analyzing Landsat Image Metadata with the Spatial DataFrame

A guest post by Gregory Brunner

A few weeks ago, Esri released an update to the ArcGIS API for Python. The newest release includes:

Hopefully, you can tell that the new functionality in the API that I am most excited about is the spatial dataframe! The spatial dataframe extends the pandas dataframe by adding geometry, spatial reference, and other spatial components to the dataframe. In adding the spatial dataframe to the API, ArcGIS users can now read feature classes, feature services, and image services directly into a dataframe. Once in a spatial dataframe, users can perform fast statistical and spatial analysis on the data, update existing feature services, and convert the dataframe to a feature class or shapefile. These are just a few examples of how you can use the spatial dataframe.

What really interests me is how this can be used with an ArcGIS image service. Can I use the spatial dataframe to extract image footprints from an image service? Can I use it to perform statistical analysis image footprints over a specific part of the world?

The answer to both of these questions is Yes! In this post, I’ll walk through how to use the API for Python to extract image service footprints from the Landsat 8 Views image service, show how to use a spatial filter to extract only footprints over New Jersey, determine the mean cloud cover and most recent acquisition date of the images, and share those image footprints as a feature service. If you have ever been interested in doing any of these, keep reading!

Getting Started

When using the ArcGIS API for Python, we first need to import the GIS.

import arcgis
from arcgis.gis import GIS
from IPython.display import display

# create a Web GIS object
gis = GIS("https://www.arcgis.com", "gregb")

We are going to look at Landsat footprints over the United States, so using gis.map we can open up our map over the USA.

#Looking over the USA
map1 = gis.map("USA", zoomlevel = 4)
display(map1)

View larger map

Using search I can find the Landsat 8 Views service in ArcGIS.com and add that layer to the map.

landsat_search = gis.content.search('"Landsat 8 Views"',outside_org=True)
map1.add_layer(landsat_search[0])

View larger map

I also need to grab the URL for this service and assign it to a variable.

url = landsat_search[0].url 
print(url)

https://landsat2.arcgis.com/arcgis/rest/services/Landsat8_Views/ImageServer

I am going to read the Landsat 8 Views service into a spatial dataframe. In order to do so, I will import SpatialDataFrame and ImageryLayer from arcgis.

from arcgis import SpatialDataFrame
from arcgis.raster._layer import ImageryLayer

Using ImageryLayer I will create the image service object.

imserv = ImageryLayer(url=url)

Within the image service object is a spatial dataframe. I can access it as follows.

image_service_sdf = imserv.query().df

I did not submit any query parameters here. I will do that shortly to show how you can specify query parameters such as a spatial filter.

The dataframe is essentially the attribute table.

image_service_sdf.head()

Let’s see how many footprints are in the dataframe.

print("There are " + str(len(image_service_sdf)) + " Landsat image service footprints in this dataframe.")

There are 1000 Landsat image service footprints in this dataframe.

There are 1000 Landsat image footprints in this dataframe. There are really hundreds of thousands, if not millions of images in this service. Only 1000 are returned because the service has a Max Record Count that is set to 1000.

Applying Query Parameters

Now that I know how to get the image service table as a spatial dataframe, I will apply a spatial filter so that I only retrieve image service footprints over New Jersey. I will also specify a where clause so that I only retrieve Primary footprints, meaning that I will exclude any Overview footprints from the dataframe.

In order to specify my extent, I am going to use arcpy and the Describe method to read a feature class that holds the geometry for New Jersey.

import arcpy

fc = r'C:\PROJECTS\STATE_OF_THE_DATA\DATA\usa_states\usa.gdb\nj'

grid_desc = arcpy.Describe(fc)
grid_sr = grid_desc.spatialReference
grid_extent = grid_desc.extent

In order to use this extent as a spatial filter on my image service query, I need to import filters from the ArcGIS API for Python.

from arcgis.geometry import filters

To create my filter, I need to pass the filter the geometry, the spatial reference, and the spatial relationship that I want the extent and the image service to have. In this case, I want to return footprints that intersect with New Jersey.

geometry = grid_extent
sr = grid_sr
sp_rel = "esriSpatialRelIntersects"
sp_filter = filters._filter(geometry, sr, sp_rel)

I also only want to return Primary image service footprints, meaning that I want to exclude the image service Overviews. I can do this by querying the Category field in the image service attribute table for footprints that have a value of 1.

search_field = 'Category'
search_val = 1
where_clause = search_field + '=' + str(search_val)

Now, when I query the image service and return the dataframe, I will also submit the where_clause and sp_filter. This should return only Primary image service footprints over New Jersey.

image_service_sdf = imserv.query(where=where_clause,geometry_filter=sp_filter,return_geometry=True).df
image_service_sdf.head()

In order to verify this, I can check and see how many footprints were returned.

print("There are " + str(len(image_service_sdf)) + " Landsat image service footprints over New Jersey.")

There are 535 Landsat image service footprints over New Jersey.

Viewing the Footprints in ArcGIS Online

Before I can convert the image service footprints to a feature class, I need to convert the AcquisitionDate from Unix time to a pandas datetime object; otherwise, the API for Python will throw an error. This can be done very easily with the pandas to_datetime method.

import pandas as pd
image_service_sdf['AcquisitionDate'] = pd.to_datetime(image_service_sdf['AcquisitionDate'] /1000,unit='s')

Now, I can verify the spatial filter worked by viewing the footprints to ArcGIS Online or viewing them in ArcMap. I will add them to ArcGIS Online using import_data.

item = gis.content.import_data(df=image_service_sdf)

I will verify that the item now exists in portal as a Feature Layer.

item

I will also load the footprints onto the map to visualize them.

map2 = gis.map("New Jersey", zoomlevel = 6)
display(map2)
map2.add_layer(item, {"renderer":"None"})

View larger map

Exporting Footprints to a Feature Class

The spatial dataframe can very easily be converted into a feature class using the to_featureclass method. I will save the image service footprints to a feature class named nj_landsat_footprints.

import os
footprint_fc =  r'C:\PROJECTS\STATE_OF_THE_DATA\DATA\usa_states\usa.gdb\nj_landsat_footprints'
image_service_sdf.to_featureclass(os.path.dirname(footprint_fc),
                                os.path.basename(footprint_fc),
                                overwrite='overwrite')

Performing Statistical Calculations on the Landsat Metadata

My motivation for doing this analysis is actually to demonstrate how easily I can use the spatial dataframe to perform statistical summaries of the imagery metadata. Having created a dataframe of only imagery metadata over New Jersey, I can now do statistical calculations that give us additional insight into our imagery data.

Collection Date

I am interested in how recently imagery over New Jersey was collected. I can use max to find the most recent image dates from my dataframe.

import datetime
most_recent_image_date = image_service_sdf['AcquisitionDate'].max()
print(most_recent_image_date)
#datetime.datetime.utcfromtimestamp(most_recent_image_date/1000).strftime('%Y-%m-%dT%H:%M:%SZ')

2017-07-14 15:40:05.235000

Similarly, I can use min to find the oldest image date in the dataframe.

mean_image_date = image_service_sdf['AcquisitionDate'].min()
print(mean_image_date)
#datetime.datetime.utcfromtimestamp(mean_image_date/1000).strftime('%Y-%m-%dT%H:%M:%SZ')

2013-08-20 15:41:30.238000

Cloud Cover

I can also analyze cloud cover in the Landsat scenes. What is the mean cloud cover of the scenes over New Jersey? How many scenes are there with zero cloud cover? These are questions that we can quickly answer using the spatial dataframe.

mean_cloud_cover = image_service_sdf['CloudCover'].mean()
print(mean_cloud_cover)

0.440778691589

The mean Cloud Cover is .4407.

I can also use the dataframe to find out how many scenes have fewer than 10% Cloud Cover.

pd.value_counts(image_service_sdf['CloudCover'].values<0.1, sort=False)

There are 131 scenes of the 535 (about 25% of scenes) that have less than 10% Cloud Cover.

Hopefully you can see the usefulness of using a spatial dataframe with image services. I’m just scratching the surface here. If you have any questions, let me know!

Published Web App

Last week Missouri got a lot of rain, which caused flooding in the Meramec river valley near St. Louis. I received a link to an ArcGIS Image Service that Surdex had created and provided that showed the flooding as of Tuesday 5/2. They flew this imagery as a service to the government and first responders, hoping to help those affected quicker and more efficiently. Here’s a sample:

I was really interested in this imagery, but it lacked some context – I couldn’t tell how the flooding depicted in the imagery compared to normal conditions. I knew ArcGIS Web AppBuilder had a swipe widget though, so I quickly:

  1. Added the Image Service (provided by Surdex) to a web map, and chose the ArcGIS Imagery basemap and an initial extent
  2. Created a Web App from that web map using Web AppBuilder.
  3. In Web AppBuilder, enabled the swipe tool to be enabled by default.
  4. Shared the web map and app publicly.

This app provided a really compelling, interesting view of the data:


(full app)

The app web viral on Twitter and Facebook – gaining 20,000 views over a three-day period, and even got published on the St. Louis Post Dispatch.

So I’m now an official, published map author now, right? 🙂

A few links:

St. Louis Lambert Airport Passenger Data

Our city’s airport posts data on how many passengers travel through the airport each month and year. It’s posted as a PDF but the data is hard to use in that format. I grabbed the data and put it into a better format:

Why not convert to a graph too?

I’ve posted this data on GitHub. It tells an interesting story of how the early 2000s was tough years for Lambert, but recent years the demand has ticked up. What other insights can we glean from better open data?

Open Data in University City

I live in University City, MO, which has no crime. Don’t believe me? This map says so:

crime-reports-u-city

Open Government

The benefits of open government data – crime data being just one aspect of that – is well documented. It allows you and I to take the data the process it, innovate on it, be informed based on it, and keep our government officials accountable.

It is easier now more than ever to have a basic open data program at any government level. Your government can use a product like ArcGIS Open Data or Socrata which are both fancy and have a lot of nice open government features. Or for cheaper, a government representative can create a free GitHub account and upload the data. Or even cheaper AND easier – simply post CSV files on an existing website. The key is that the data be machine readable and in an open format.

University City was not doing that – crime data, before a few months ago, was being released in non-machine-readable PDFs. Thus it’s not able to be placed on crime maps, or analyzed by experts, or any of that. And so we end up with the hole in the crime map over our town.

“Let’s fix this!” I said. Little did I know it’d be a years-long adventure.

The initial contact

I originally reached out to my city council representative, mentioning the benefits of open data and why our city should be doing more to provide open data.

10/27/2013
To: City Council Rep
From: Gavin

I was just wondering if you know of any current or future plans of U City to open our city data? I recently wrote a post on open data – http://www.gavinr.com/2013/10/16/share-data-github/ – and would love to help U City in getting our public data more available in open formats.

I got a quick response saying my message was received, but after a few back-and-forth emails, nothing ever happened.

Go to the Top!

A couple years later, I emailed the police chief directly:

3/26/2015
To: Colonel Adams (U City Police Chief)
CC: City Council Rep
From: Gavin

I am writing to inquire why the crime stats posted on http://www.ucitymo.org/index.aspx?nid=482 are posted in a hard-to-process PDF format, instead of a more open CSV (Excel) format. I think the city at large would benefit greatly from having that data more easily process-able. Is there any way your team could provide this?

Response from this: nothing. Crickets. I took this to mean neither my city council rep nor the police chief has any interest in helping me get this done.

Finally

In 2016, University City adopted “MyGov” as their website CMS. One portion of this solution allows citizens to make a request digitally to the city. One of those options is “Police Issues – Crime Statistics.” Interesting. So I fill out the form, using the same request as above and wait. No response.

XLSX files on University City homepage

The new files, unceremoniously relegated to the bottom of the crime stats page

After three months of waiting, I forwarded the email to my city council rep, who forwarded to the city manager, who forwarded to the police department, who finally acknowledged my request and posted the crime data in XLSX format (only this year). I had to further ask for previous years to be posted.

WHEW. Finally after almost three years of discussing this, University City has made a small step toward open data. I learned through this process that when making requests to government officials, you sometimes have to be extremely specific, because general ideas or initiatives will probably get dropped. You also need to persevere – doing the right thing is not always the easiest thing (although it’s quite easy now – uploading an XLSX file? c’mon), so you’ll probably have to ask multiple times, and multiple ways, to get them to do what’s right.

What’s Next?

The good news: my city is posting the data in a decent format. Bad news: it’s not the best layout of the data. In a later post I’ll show how I wrote some scripts to process that and get it into databases and maps.

Extending Widgets in ArcGIS Web AppBuilder

I’ve had many requests and ideas for Web AppBuilder widgets over the past few years, and many times they are ideas that are building on an existing widget. For example, a query widget that is custom to a particular workflow, an edit widget that enforces certain geometry movement restrictions, or maybe adding a few more options to the built-in print widget.

Most developers will copy the original out-of-the-box (OOTB) widget, rename it, and then start development from there. Which is great because it gets you moving quickly. But is this really the best way to do things?

There are many problems that you’ll eventually encounter when doing this copy/paste dance, and they all stem from the fact that there is no separation of your custom code from the original, OOTB widget code. Some potential future problems include:

  • Upgrading your widget to a new version of the OOTB widget is almost impossible
  • It’s hard to identify if bugs you encounter are stemming from your custom code or the original widget code.
  • Future developers trying to understand how your widget works have a much larger codebase to learn and understand (your custom code plus the original widget code).

Luckily there is a good way of “extending” modules in dojo/declare. Most widgets use this feature to extend from jimu/BaseWidget, but you can use it to extend from an existing widget too!

So how would you do this? (I’ll use this repo as an example)

First create a new, blank widget. You can do that by either creating a manifest.json and Widget.js file manually, or use the convenient yeoman generator that will create the files automatically.

Next, modify your widget to extend from an existing widget instead of extending from BaseWidget. You have to include the “parent widget” in the page as well as referencing it in your delcare statement. See an example here.

Next, include the translation and style files from the original widget. Don’t forget to set this.nls.

Finally, override functions in the parent widget as needed to get your desired functionality. In this example, I overrode the _createEditor function, but this will be different depending on what you want to do with your widget.

So happy from my clean widget code!

Extending your widget is a little more work, but you end up with much cleaner code. I hope you give it a try!

Have you done this or something similar? Let me know in the comments, or contribute to the example GitHub repository. Want more Web AppBuilder development tips and tricks? Catch my talk at Esri Dev Summit 2017!

Older Posts

Gavin Rehkemper

JavaScript, WordPress, and GeoDev

© 2017 Gavin Rehkemper | JavaScript, WordPress, and GeoDev