Skip navigation
All Places > Applications Prototype Lab > Blog > Author: david_johnson-esristaff


In June of 2017 we began another collaboration with Dr. Camilo Mora of the University of Hawaii, Department of Geography. This came on the heels of our previous project with Dr. Mora to develop a web mapping application to display his team's research on climate change and deadly heatwaves. For their next project they had expanded their research to include multiple cumulative hazards to human health and well-being resulting from climate change.  These hazards include increased fires, fresh water scarcity, deforestation, and several others. Their new research was recently published in the journal Nature Climate Change.  Several news outlets published stories on their findings, including these from The New York Times, Le Monde, and Science et Avenir. For our part, the Applications Prototype Lab developed an interactive web mapping application to display their results. To view the application, click on the following image. To learn how to use the application, and about the research behind it, click on the links for "Help" and "Learn More" at the top of the application.


Cumulative Exposure to Climate Change


In this post I'll share some of the technical details that went into the building of this application.  


The Source Data


For each year of the study, 1956 - 2095, the research team constructed a series of global data sets for 11 climate-related hazards to human health and well-being. From those data sets they built a global cumulative hazards index for each year of the study. For information about their methods to generate these data sets, refer to their published article in Nature Climate Change. Each data set contains the simulated (historical) or projected (future) change in intensity of a particular hazard relative to a baseline year of 1955. For the years 2006 - 2095, hazards were projected under three different scenarios of greenhouse gas (GHG) emissions ranging from a worst-case business-as-usual scenario to a best-case scenario where humanity implements strong measures to reduce GHG emissions. In total, they produced 3828 unique global data sets of human hazards resulting from climate change.


Data Pre-processing


We received the data as CSV files which contained the hazard values on a latitude-longitude grid at a spatial resolution of 1.5 degrees. The CSV data format is useful for exchanging data between different software platforms. However, it is not a true spatial data format. So we imported the data from the CSV files into raster datasets. This is typically a two-step process where you first import the CSV files into point feature classes and then export the points to raster datasets. However, since the data values for the 11 hazards were not normalized to a common scale, we added a step to re-scale the values to a range of 0 - 1, according to the methodology of the research team, where:  

  • 0 equals no increase in hazard relative to the historical baseline value.
  • 1 equals the value at the 95th percentile or greater of increased hazard between 1955 and 2095 for the "business-as-usual" GHG emissions scenario.


With a spatial resolution of 1.5 degrees, each pixel in the output raster datasets are approximately 165 Km in width and height. This was too coarse for the web app, because the data for land-based hazards such as fire and deforestation extended quite a distance beyond the coastlines. So we added another processing step to up-sample each dataset by a factor of ten and remove the pixels from the land-based hazards raster datasets whose centers were outside of a 5 Km buffer of the coastlines.  

upsampling and clipping


We automated the entire process with Python scripts, using geoprocessing tools to convert the source data from CSV to raster dataset, build the coastal buffer, and up-sample and clip the land raster datasets.  To re-scale the data values, we used mathematical expressions. At the end of these efforts we had two collections of raster datasets - one for the 11 hazards indexes, and another for the cumulative hazards index.


Data Publishing


We built two mosaic datasets to organize and catalog each collection of raster datasets. From each mosaic dataset we published an image service to provide the web application with endpoints through which it could access the data. On the application, the map overlay layer is powered by the image service for the cumulative hazards index data. This layer is displayed in red with varying levels of transparency to indicate the level of cumulative hazards at any location. To support this type of rendering, we added a custom processing template to the image service's source mosaic dataset. The processing template uses the Stretch function to dynamically re-scale the floating-point data values in the source raster datasets to 8-bit integers, and the Attribute Table function to provide the color and transparency values of the exported image on a per-pixel basis.


The Animations


We built short video animations of the change in cumulative hazards over time using the Time and Animation Toolbars in ArcGIS Pro. You can access those animations from the application by clicking on the "Animations" link at the top of the application window. We used the cumulative hazards index image service as the data source of the animation. This service is time-aware, enabling us to define a timeline for the animations. Using the capabilities in the Animations Toolbar, we defined properties such as the time-step interval and duration, total duration, output format and resolution, and the various overlays such as the legend, watermarks, and dynamic text to display the year. We filtered the data in the image service by GHG emissions scenario using definition queries to create three separate animations of the change in cumulative hazards over time.


The Web Application


We built the web application using the ArcGIS API for JavaScript. To render the cumulative hazards map layer, the application requests the data from the image service in the LERC format.  This format enables the application to get the color and transparency values for each pixel from the attribute table to build a client-side renderer for displaying the data. The chart that appears when you click on the map was built with the Dojo charting library. This chart is powered by the image service with the 11 individual human hazards index data. To access the hazards data, the web application uses the Identify function to get the values for each of the 11 hazards at the click location with a single web request to the service. 


In Summary


Building this application gave us the opportunity to leverage many capabilities in the ArcGIS platform that are well suited for scientific analysis and display. If you are inspired to build similar applications, then I hope this post provides you with some useful hints. If you have any technical questions, add them into the comments and I'll try to answer them. I hope this application helps to extend the reach of this important research as humanity seeks to understand the current and projected future impacts of climate change.

CityEngine Station-Hand model


Every now and then a really unique and out-of-the-box idea comes our way that expands our conceptions about the possible applications of the ArcGIS platform. This was one of those ideas. Could GIS be used to map the human body? More specifically, could we use CityEngine to visualize the progress of physical therapy for our friend and Esri colleague Pat Dolan of the Solutions team? Pat was eager to try this out, and he provided a table of measurements taken by his physical therapist to track his ability to grip and extend his fingers over time. With the help of the CityEngine team, we developed a 3D model of a hand, and used CityEngine rules to apply Pat's hand measurements to the model. We thought it would be fun to show a train station in a city that would magically transform into a hand. Our hand model is not quite anatomically correct, but it has all the digits and they are moveable!


Click the image above to view a short video of this project. Pat and I showed this application, and others, at the 2017 Esri Health and Human Services GIS Conference in Redlands. Click here to view a video of that presentation.


Among the best resources for learning the ArcGIS API for Python are the sample notebooks at the developers website. A new sample notebook is now available that demonstrates how to perform a network analysis to find the best locations for new health clinics for amyotrophic lateral sclerosis (ALS) patients in California. To access the sample, click on the image at the top of this post.


I originally developed this notebook for a presentation that my colleague Pat Dolan and I gave at the Esri Health and Human Services GIS Users conference in Redlands, California in October. Although network analysis is available in many of Esri's offerings, we chose the Jupyter Notebook, an open-sourced browser-based coding environment, to show the attendees how they could document and share research methodology and results using the ArcGIS API for Python.  This sample notebook provides a brief introduction to network analysis and walks you through our methodology for siting new clinics, including accessing the analysis data, configuring and performing analyses, and displaying the results in maps. 

One of the great things about working in the Lab is you get to experiment with the new goodies from our core software developers before they are released.  When I heard that version 1.2 of the ArcGIS API for Python would include a new module for raster functions, I could not wait to give it a try.  Now that v.1.2 of the API is released, I can finally show you a Jupyter Notebook I built which has an example of a weighted overlay analysis implemented with raster functions.   The following is a non-interactive version of that notebook which I exported to HTML.  I hope it will give you some ideas for how you could use the ArcGIS API for Python to perform your own raster analysis.



Finding Natural and Accessible Areas in the State of Washington, USA

The weighted overlay is a standard GIS analysis technique for site-suitability and travel cost studies. This notebook leverages the new "arcgis.raster.functions" module in the ArcGIS API for Python 1.2 to demonstrate an example of a weighted overlay analysis.  This example attempts to identify areas in the State of Washington that are "natural" while also being easy to travel within based on the following criteria:

  • elevation (lower is better)
  • steepness of the terrain (flatter is better)
  • degree of human alteration of the landscape (less is better)

The input data for this analysis includes a DEM (Digital Elevation Model), and a dataset showing the degree of human modification to the landscape.

In general, weighted overlay analysis can be divided into three steps:

  1. Normalization: The pixels in the input raster datasets are reclassified to a common scale of numeric values based on their suitability according to the analysis criteria.
  2. Weighting: The normalized datasets are assigned a percent influence based on their importance to the final result by multiplying them by values ranging from 0.0 - 1.0. The sum of the values must equal 1.0.
  3. Summation: The sum of the weighted datasets is calculated to produce a final analysis result.


We'll begin by connecting to the GIS and accessing the data for the analysis.

Connect to the GIS

In [1]:
# import GIS from the arcgis.gis module
from arcgis.gis import GIS

# Connect to the GIS.
   web_gis = GIS("", 'djohnsonRA')   
   print("Successfully connected to {0}".format(
Enter password:········ 
Successfully connected to ArcGIS Enterprise A

Search the GIS for the input data for the analysis

Human Modified Index

In [2]:
# Search for the Human Modified Index imagery layer item by title
item_hmi ='title:Human Modified Index', 'Imagery Layer')[0]
Human Modified Index 
A measure of the degree of human modification, the index ranges from 0.0 for a virgin landscape condition to 1.0, for the most heavily modified areas.Imagery Layer by djohnsonRA 
Last Modified: July 06, 2017 
0 comments, 2 views


In [3]:
# Search for the DEM imagery layer item by title
item_dem ='title:USGS NED 30m', 'Imagery Layer')[0]
The National Elevation Dataset (NED) is the primary elevation data product of the USGS. This version was resampled to 30m from source data at 1/3 arc-second resolution and projected to an Albers Equal Area coordinate system.Imagery Layer by djohnsonRA 
Last Modified: July 06, 2017 
0 comments, 8 views

Study area boundary and extent

In [4]:
# Search for the Ventura County feature layer item by title
item_studyarea ='title:State of Washington, USA',
                                        'Feature Layer')[0]
State of Washington, USA 
State of Washington, USAFeature Layer Collection by djohnsonRA 
Last Modified: July 07, 2017 
0 comments, 2 views
In [5]:
# Get a reference to the feature layer from the portal item
lyr_studyarea = item_studyarea.layers[0]

Get the coordinate geometry of the study area

In [6]:
# Query the study area layer to get the boundary feature
query_studyarea = lyr_studyarea.query(where='1=1')
# Get the coordinate geometry of the study area.
# The geometry will be used to extract the Elevation and Human Modified Index data.
geom_studyarea = query_studyarea.features[0].geometry
# Set the spatial reference of the geometry.
geom_studyarea['spatialReference'] = query_studyarea.spatial_reference

Get the extent of the study area

In [7]:
# Import the geocode function
from arcgis.geocoding import geocode
# Use the geocode function to get the location/address of the study area
geocode_studyarea = geocode('State of Washington, USA',
out_sr= query_studyarea.spatial_reference)
In [8]:
# Get the geographic extent of the study area
# This extent will be used when displaying the Elevation, Human Modified Index,
# and final result data.
extent_studyarea = geocode_studyarea[0]['extent']
{'xmax': -1451059.3770040546,  
'xmin': -2009182.5321227335, 
'ymax': 1482366.818700374, 
'ymin': 736262.260048952}

Display the analysis data

Human Modified Index

In [9]:
# Get a reference to the imagery layer from the portal item
lyr_hmi = item_hmi.layers[0]
# Set the layer extent to geographic extent of study area and display the data.
lyr_hmi.extent = extent_studyarealyr_hmi


In [10]:
# Get a reference to the imagery layer from the portal item
lyr_dem = item_dem.layers[0]
# Set the layer extent to the geographic extent of study area and display the data.
lyr_dem.extent = extent_studyarealyr_dem

Slope (derived from elevation via the Slope raster function)

In [11]:
# Import the raster functions from the ArcGIS API for Python (new to version 1.2!)
from arcgis.raster.functions import *
In [12]:
# Derive a slope layer from the DEM layer using the slope function
lyr_slope = slope(dem=lyr_dem,slope_type='DEGREE', z_factor=1)
# Use the stretch function to enhance the display of the slope layer.
lyr_slope_stretch = stretch(raster=lyr_slope, stretch_type='StdDev', dra='true')
# Display the stretched slope layer within the extent of the study area.
lyr_slope_stretch.extent= extent_studyarealyr_slope_stretch

Extract the data within the study area geometry

Use the Clip raster function to extract the analysis data from within the study area geometry

Human Modified Index

In [13]:
# Extract the Human Modified Index data from within the study area geometry
hmi_clipped = clip(raster=lyr_hmi, geometry=geom_studyarea)


In [14]:
# Extract the Elevation data from within the study area geometry
elev_clipped = clip(raster=lyr_dem, geometry=geom_studyarea)


In [15]:
# Extract the Slope data from within the study area geometry
slope_clipped = clip(raster=lyr_slope, geometry=geom_studyarea)
# Apply the Stretch function to enhance the display of the slope_clipped layer.
slope_clipped_stretch = stretch(raster=slope_clipped, stretch_type='StdDev',

Perform the analysis

Step 1: Normalization

Use the Remap function to normalize each set of input data to a common scale of 1 - 9, where 1 = least suitable and 9 = most suitable.

In [16]:
# Create a colormap to display the analysis results with 9 colors ranging 
# from red to yellow to green.
clrmap=  [[1, 230, 0, 0], [2, 242, 85, 0], [3, 250, 142, 0], [4, 255, 195, 0],
         [5, 255, 255, 0], [6, 197, 219, 0], [7, 139, 181, 0], [8, 86, 148, 0],
9, 38, 115, 0]]
In [17]:
# Normalize the elevation data
elev_normalized = remap(raster=elev_clipped,
                        input_ranges=[0,490, 490,980, 980,1470, 1470,1960, 1960,2450,
                                      2450,2940, 2940,3430, 3430,3700, 3920,4100],
                       output_values=[9,8,7,6,5,4,3,2,1], astype='U8')

# Display color-mapped image of the reclassified elevation data
colormap(elev_normalized, colormap=clrmap) 
In [18]:
# Normalize the slope data
slope_normalized = remap(raster=slope_clipped,                         
                        input_ranges=[0,1, 1,2, 2,3, 3,5, 5,7, 7,9, 9,12, 12,15,
                        output_values=[9,8,7,6,5,4,3,2,1],  astype='U8') 

# Display a color-mapped image of the reclassified slope data
colormap(slope_normalized, colormap=clrmap)
In [19]:
# Normalize the Human Modified Index data
hmi_normalized = remap(raster=hmi_clipped,                 
                      input_ranges=[0.0,0.1, 0.1,0.2, 0.2,0.3, 0.3,0.4, 0.4,0.5,
0.5,0.6, 0.6,0.7, 0.7,0.8, 0.8,1.1],                 
                      output_values=[9,8,7,6,5,4,3,2,1],  astype='U8')

# Display a color-mapped image of the reclassified HMI data
colormap(hmi_normalized, colormap=clrmap)

Step 2: Weighting

Use the overloaded multiplication operator * to assign a weight to each normalized dataset based on their relative importance to the final result.

In [20]:
# Apply weights to the normalized data using the overloaded multiplication 
# operator "*".
# - Human Modified Index: 60%
# - Slope: 25%
# - Elevation: 15%
hmi_weighted = hmi_normalized * 0.6
slope_weighted = slope_normalized * 0.25
elev_weighted = elev_normalized * 0.15

Step 3: Summation

Add the weighted datasets together to produce a final analysis result.

In [21]:
# Calculate the sum of the weighted datasets using the overloaded addition 
# operator "+".
result_dynamic = colormap(hmi_weighted + slope_weighted + elev_weighted,
                          colormap=clrmap, astype='U8')

The same analysis can also be performed in a single operation

In [22]:
result_dynamic_one_op = colormap(    
# Human modified index layer       
      0.60 * remap(raster=clip(raster=lyr_hmi, geometry=geom_studyarea),
input_ranges=[0.0,0.1, 0.1,0.2, 0.2,0.3, 0.3,0.4, 0.4,0.5,
0.5,0.6, 0.6,0.7, 0.7,0.8, 0.8,1.1],                    
      # Slope layer       
      0.25 * remap(raster=clip(raster=lyr_slope, geometry=geom_studyarea),
input_ranges=[0,1, 1,2, 2,3, 3,5, 5,7, 7,9, 9,12, 12,15,
      # Elevation layer       
      0.15 * remap(raster=clip(raster=lyr_dem, geometry=geom_studyarea),
         input_ranges=[-90,250, 250,500, 500,750, 750,1000, 1000,1500,
                       1500,2000, 2000,2500, 2500,3000, 3000,5000],                   
   colormap=clrmap,  astype='U8')

Generate a persistent analysis result via distributed server based raster processing.

Portal for ArcGIS has been enhanced with the ability to perform distributed server based processing on imagery and raster data. This technology enables you to boost the performance of raster processing by processing data in a distributed fashion, even at full resolution and full extent.

You can use the processing capabilities of ArcGIS Pro to define the processing to be applied to raster data and perform processing in a distributed fashion using their on premise portal. The results of this processing can be accessed in the form of a web imagery layer that is hosted in their ArcGIS Organization.

For more information, see Raster analysis on Portal for ArcGIS

In [23]:
# Does the GIS support raster analytics?
import arcgis
In [24]:
# The .save() function invokes generate_raster from the
# module to run the analysis on a GIS server at the source resolution of the
# input datasets
and store the result as a persistent web imagery layer in the GIS.
result_persistent ="NaturalAndAccessible_WashingtonState")
Analysis Image Service generated from GenerateRasterImagery Layer by djohnsonRA 
Last Modified: July 07, 2017 
0 comments, 0 views
In [25]:
# Display the persistent result
lyr_result_persistent = result_persistent.layers[0]
lyr_result_persistent.extent = extent_studyarea
Data Credits:
A measure of the degree of human modification, the index ranges from 0.0 for a virgin landscape condition to 1.0 for the most heavily modified areas. The average value for the United States is 0.375. The data used to produce these values should be both more current and more detailed than the NLCD used for generating the cores. Emphasis was given to attempting to map in particular, energy related development. Theobald, DM (2013) A general model to quantify ecological integrity for landscape assessment and US Application. Landscape Ecol (2013) 28:1859-1874 doi: 10.1007/s10980-013-9941-6
USGS NED 30m:  
Data available from the U.S. Geological Survey. See USGS Visual Identity System Guidance for further details. Questions concerning the use or redistribution of USGS data should be directed to: or 1-888-ASK-USGS (1-888-275-8747). NASA Land Processes Distributed Active Archive Center (LP DAAC) Products Acknowledgement: These data are distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at USGS/EROS, Sioux Falls, SD.
State of Washington: Esri Data & Maps

Annual number of lethal heat days if global carbon emissions are not reduced

In the Spring of last year Dr. Camilo Mora of the University of Hawaii Manoa contacted our team.  He wanted to know if we would be interested in developing an interactive web map to display the results of a research project he was leading into the effect of rising global temperatures on climatic conditions resulting in human mortality due to heat stress.  We were glad to hear from Dr. Mora again.  The year before we had developed a web map to display the results of his research into how climate change could affect the length of plant growing seasons.  For this new study, Dr. Mora’s team analyzed 1763 documented lethal heat events and identified a global threshold beyond which daily mean surface air temperature and relative humidity became deadly.  Using Earth Systems Models to predict daily values for temperature and humidity across the globe up to the year 2100, they estimated the likely number of lethal heat days annually under low, moderate, and high carbon emissions scenarios.


Although his new research project was still in it's early stages, we found the initial results to be very compelling and we agreed to move forward with the project.  Using preliminary data from their research, we explored some ideas for how to present the data and developed a couple of prototype applications.  Several months later, we heard from Dr. Mora again. His team had completed their research, and he was ready to share his finalized data with us and to collaborate on the design of the final application.   The time-frame was short.  Dr. Mora and his team were writing the final drafts of a paper for publication in the journal Nature Climate Change.   So we rolled up our sleeves, reviewed our initial prototypes, explored the finalized data, and then got straight to work. 


The application Heatwaves: Number of deadly heat days leverages the robust capabilities of the ArcGIS platform to distill complex scientific data into intuitive maps that enable users to interact with and understand the data.  This was an interesting development project, not only for it's subject matter, but also on a technical level.  So we thought it would be worthwhile to share some details about how we built the application.


Heatwaves: Number of deadly heat days


At the front-end of the application is a web map developed using jQuery, Bootstrap, and the ArcGIS API for Javascript.  The map contains an image service layer which displays the number of lethal heat days at any location over land using a color ramp with a range from white to yellow to red to black representing 0 - 365 days respectively for any year from 1950 - 2100.  You can select the year from a slider control or from a drop down list.  The data on the annual number of lethal heat days for the years 1950 – 2005 are based on historic climate records. 

RCP List

The data for the years 2006 - 2100 are based on the average of 20 Earth System Models developed for the Coupled Model Intercomparison Project Phase 5, under low, moderate, and high (i.e. "business as usual") carbon emissions scenarios (i.e. Representative Concentration Pathways, RCPs 2.6, 4.5, and 8.5 respectively).  By selecting from a drop-down list of RCPs you can view the modeled results for the different carbon emissions scenarios.

Heatwaves App Identify WindowWhen you click on a location over land, a window appears with a line chart and a scatter plot that reveal further insights into the study results for that location.  The line chart displays the trend in the annual number of lethal heat days at the location for each year of the study period.  The scatter plot displays the temperature and humidity for each day of the selected year over a curve which represents the lethal heat threshold.

Now let's take a look at some of the deeper technical details of this application.  On the back-end of the application are two web services that deliver the data from the study results to the web application for display.  These services are hosted on the Esri Applications Prototype Lab's GIS Portal. 

An image service provides the web application with the data for the annual number of lethal heat days for each year of the study period.  The data source of the service is a mosaic dataset that defines a single point of access to a collection of single-band raster datasets of global extent.  Each raster dataset contains the number of lethal heat days across the globe for a given year.  For the historical period 1950–2005, the data for each year are stored in a single raster dataset.  For the future period 2006–2100, the data for each year are stored in 3 raster datasets – one for each of the carbon emissions scenarios. 

The image service has two roles: 1). to provide the images showing the annual number of lethal heat days for display in the web map  2). to provide the data for the graph of the trend in time of the annual number of lethal heat days.  To generate the images for the map layer, the mosaic dataset applies a raster function chain that dynamically clips the source raster datasets to the coastlines and applies a color ramp to convert the single-band source data into three-band color RGB output images. To provide the data for the trend graph, the service delivers the pixel values at a given location from each of the historic rasters and from the future rasters for the selected carbon emissions scenario.

A geoprocessing service provides the data for the chart that plots the temperature and relative humidity for each day of a given selected year.  The source data for this service are a collection of 36 NetCDF files that contain the daily values for temperature and relative humidity for the study period and for each carbon emissions scenario.  Each file contains data for a twenty year period for either temperature or relative humidity for the historic period, or for one of the three carbon emissions scenarios.  In total, the files use 17 GB of storage and contain 12,570,624 unique points of data.  To build this service, we started by developing a Python script with input parameters for the selected year, the selected carbon emissions scenario, and the coordinates of the location where the user clicked on the map.  The script obtains the requested data from those files in four steps:

  1. The NetCDF files containing the relevant temperature and humidity data are identified from the first two input parameters.  
  2. In-memory tables are created from the files using the Make NetCDF Table View geoprocessing tool.
  3. Queries are crafted to obtain the temperature and humidity values from the tables for each day of the selected year at the specified location. 
  4. The results of the queries are sorted by day of the year and returned the to client application.

The python script was then wrapped into a Python script tool and published as a geoprocessing service.

The application also includes links three video animations showing the increase in lethal heat days over time for each of the carbon emissions scenarios.  These videos were created using the animation tools in ArcGIS Pro.  The videos representing rcp2.6, rcp4.5, and rcp8.5 can be viewed here, here, and here.  Links to the videos and the source code of the application are also available from the application when you click the "About" button at the top right corner.

In conclusion, we'd like to thank Dr. Mora and his team for their very important research and for the opportunity to contribute in our way towards helping to extend the reach of their findings.  We enjoyed working with Dr. Mora and hope to collaborate with him on his future projects.


In the news

Deadly heat waves becoming more common due to climate change | CNN

Deadly Heat Waves Could Threaten 3 in 4 People by 2100 | HUFFPOST

Half of World Could See Deadly Heat Waves By 2100 | Climate Central

Study shows deadly heat waves are becoming more frequent | Chicago Tribune

A third of the world now faces deadly heatwaves as result of climate change | The Guardian

By 2100, Deadly Heat May Threaten Majority of Humankind | National Geographic

Deadly heatwaves could affect 74 percent of the world's population | University of Hawaii News

Deadly Heat Waves Could Endanger 74% of Mankind by 2100, Study says | inside climate news

Billions to Face 'Deadly Threshold' of Heat Extremes by 2100, Study Finds | EcoWatch

Killer Heat Waves Will Threaten Majority of Humankind by Century's End | Alternet


In my previous post about the Combined Field Statistics addin for ArcMap I mentioned that it was inspired from another addin toolkit I am currently developing which analyzes and displays rates of disease occurrence in a map.  In this post I want to share with you another addin project which was borne out of the disease mapping project.  This one is called Diverging Color Ramp.  This addin enables you to intelligently apply a dichromatic color ramp to feature layers containing points, lines, or polygons in ArcMap 10.4 and above. It is somewhat similar in capability to the smart mapping tools in ArcGIS Online and is useful if you need to display data with values that diverge from a specific threshold value which you define.  The screenshot at the top of this post shows a map layer of population density where the colors diverge from an arbitrarily defined threshold value of 40,000/square mile.   Other examples of threshold values might include 1.0 for ratios, 0.0 for values which represent gain or loss from an average value, or any value which represents a significant threshold of interest in your data.

This tool gives you a great deal of control over how a feature layer is rendered.  It provides a pair of default color ramps, one for for values above the threshold and another for values below the threshold, as well as a separate threshold color symbol.  If you don't like the colors you can change them to anything you like by accessing them from the Style Manager in ArcMap.  You can also choose the classification method and number of classes, and even choose the number classes which are generated for values above and below the threshold.  This addin also has a "secret" trick up it's sleeve.  When it generates the symbols for each class in the map layer's renderer, it actually makes a clone of the symbol for the first class in the current renderer to use as a template for the rest of the symbols.  After that is simply sets the background color of each symbols  to colors pulled from the color ramps.  This makes it easy to quickly apply global changes to other properties in the symbols, such as outline color or fill pattern, since you only have to make the change in the first symbol and then apply it as a template to each of the symbol classes with a single button click.  This is a fun tool to use since it makes it easy to try out different ways to visualize your data.  I hope you enjoy using it.

Download links:Diverging Color Ramp (addin + documentation)Diverging Color Ramp (source code + documentation)
Combined-Field-Statistics-GUI.pngI'm currently working on a new add-in toolkit for ArcMap 10.4 to analyze and display rates of disease occurrence in a map.  This is a work-in-progress and I will provide more details about it in the coming months.  But in the meantime I want to share with you the result of a smaller project that was borne out of the disease mapping project.

One of the requirements for the disease mapping project is the ability to calculate the sums of values from age group columns in a feature layer.  Generally, the tool of choice for a task like this is the Summary Statistics geoprocessing tool in ArcToolBox.  However, I also needed the ability to calculate sums from the combined values in multiple age group columns.  For example, to calculate the total number of individuals in columns for ages 0 -4 and 5 - 9, to obtain an aggregate total for ages 0 - 9.

It occurred to me that this capability could have broader applications if the tool could also calculate the full set of statistics provided by the Summary Statistics tool.  So I decided to build a Combined-Field-Statistics-Result.pngmore generic stand-alone addin tool in parallel with the disease mapping toolkit.  The result of this effort is the Combined Field Statistics Addin.  This tool is very simple to use and includes a detailed user guide to get you started.  It also generates an output log so you can keep track of how the tool was configured to generate a particular output.  If this capability sounds useful to you give it a try and let me know what you think in the comments!

Download links:Combined Field Statistics (addin + documentation)Combined Field Statistics (source code + documentation)


The source code for Dimension Explorer, an add-in for exploring time-aware and multidimensional data, is now available for download.  Version 1.11 is a minor update from version 1.1 and is built specifically for ArcMap 10.4.  This version features tighter integration with the ArcMap Table of Contents window - as layers are added, removed or changed, the  contents of the layer selection list in the Settings window are updated accordingly.

Update - I've had reports that the layer selection list is not getting updated consistently as layers are added, removed, or had their properties updated in the map document.  Until this issue is resolved, I recommend users of ArcMap 10.3.1 and 10.4 use version 1.1 of Dimension Explorer instead.  I apologize for any inconvenience this may have caused.

Download links:

Dimension Explorer 1.11 (addin and documentation)

Dimension Explorer 1.11 (source code and documentation)

Dimension Explorer 1.1 (original blog post)

A new version of Dimension Explorer is now available for download.  Dimension Explorer 1.1 is an addin tool for exploring time-aware and multidimensional map data in ArcMap 10.3 and above.

Here is what's new in version 1.1:

    • map layers created with the Make NetCDF Raster Layer and Make NetCDF Feature Layer geoprocessing tools are now supported.


    • map layers with vertical dimension values defined as ranges (e.g. 0-10 meters, 5-50 meters, etc) are now supported.


    • export a map animation in ArcMap to a series of still images to create video animations with applications such as Windows Movie Maker and ffmpeg.


  • various bug fixes and optimizations.

Here is a video animation of the minimum Arctic sea ice extent for the years 1979 - 2014.  I created it with Windows Movie Maker using still images exported via Dimension Explorer 1.1.  The map includes a time-aware layer created with the Make NetCDF Feature Layer geoprocessing tool with data from the NOAA.

Dimension Explorer 1.1 can be downloaded here


If you are looking for data to get started with Dimension Explorer,  the NOAA Earth System Research Laboratory, Physical Sciences Division, (NOAA/ESRL PSD) has many large collections of spatial scientific data for climate and oceans in NetCDF format.  I recommend starting with their gridded climate datasets.   You can add most of their datasets to ArcMap using the Make NetCDF Raster Layer geoprocessing tool.  If you get the error "One or both dimensions have variable spacing in their coordinates", use the Make NetCDF Feature Layer geoprocessing tool instead.  If the datasets are stored as multiple files representing the data at different times, use the Mosaic Dataset to temporally aggregate the files using the NetCDF Raster Type.  Finally, if you are working with temporal data of any type, be sure to time-enable the layer in ArcMap.

DimensionExplorer4.pngUpdate - a more recent version of Dimension Explorer is now available.  Click here for more information.

Dimension Explorer, an addin tool for ArcMap, has just been released by the Esri Applications Prototype Lab!

Dimension Explorer 1.0 makes it easier to work with  time-aware and multidimensional data in ArcMap 10.3 by providing slider controls for navigation.  It works by retrieving dimensional information from a map layer to build an interactive dimensional model that can be edited and saved in the map document.  Dimension Explorer is the successor to the Timeliner addin for ArcMap, which also works in ArcMap 10.3 and can be downloaded here.

Click here to download Dimension Explorer

With the 10.3 release of ArcGIS, the mosaic dataset now supports multidimensional data in NetCDF, GRIB, and HDF format.  Dimension Explorer supports map layers based on mosaic datasets and image services which are time-aware or multidimensional, and time-aware feature classes.

Download Dimension Explorer and let us know what you think of it in the comments!

Image services are not only for serving imagery; they can also perform dynamic pixel-level analysis on multiple overlapping raster datasets using chained raster functions.  Image services deliver blazing fast performance with pre-processed source imagery, especially when it is served from a tile cache.  Dynamic image services, on the other hand, may not respond as quickly because of the additional processing demands they place on the server.  High performance is important for dynamic image services because the data is re-processed automatically each time the user pans or zooms the map.  Dynamic image services typically produce sub-second responses for simple types of analysis.  But more complex analysis may take much longer depending on the complexity of the processing chain and the number of input datasets involved.  How can we get better performance in those situations?  To answer that question I ran a series of tests to see how the following factors affect the performance of dynamic image services:


    • Map scale and source data resolution
    • Resampling method
    • Source data format and compression
    • Project on-the-fly
    • Request size

In this article I present the results of those tests along with some suggestions for maximizing performance.  My testing machine is a desktop computer running Windows 7 SP1 and ArcGIS 10.2.1 with 18GB of RAM and a quad-core Intel Xeon W3550 processor running at 3 GHZ.  The test data was stored on an otherwise empty 2 TB SATA hard drive that I defragmented and consolidated prior to testing.   The tests were configured to determine the average response times of services under various conditions.  By “response time” I mean the time it takes a service to retrieve the source data, process it, and transmit an output image.  Transmission time was minimized by running the testing application directly on the server machine.


This information is written with the intermediate to advanced GIS user in mind.  I assume the reader has a general understanding of image services, raster data and analysis, raster functions, geoprocessing, mosaic datasets, map projections, and map service caching.


Map Scale and Source Data Resolution


The pixels that are processed for analysis by dynamic image services are generally not identical to the pixels stored in the source datasets.  Instead, the source data pixels are first resampled on-the-fly to a new size based on the current scale of the map.  This formula shows the relationship between map scale and resampling size when the map units of the data are meters:


Resampled pixel size = map scale * 0.0254/96


The resampled pixel size is analogous to the “Analysis Cell Size” parameter in the Geoprocessing Framework and is sometimes referred to as the “pixel size of the request”.  As you zoom out to smaller map scales, the resampled pixel size increases until eventually the service resamples from the pixels in the pyramids.  Resampling from pyramids helps to keep the performance of the service relatively consistent over a range of map scales.



Chart 1. Performance of an image service that performs a binary overlay analysis over a range of map scales.

Performance still varies depending on map scale and typically looks similar to chart 1.    I generated these results using an application configured to simulate a single user panning the map 100 times in succession at specific map scales.  The chart shows the average time the service took to process and transmit the output images for different map scales.  This particular service was configured with a raster function template to perform a binary overlay analysis on eleven overlapping rasters in a mosaic dataset.  The pixel sizes of the source datasets ranged from 91.67 to 100 meters.  The raster function template was configured to return a binary result, where each output pixel is classified as either “suitable” or “unsuitable” based on the analysis parameters.


Take a look at the three points along the horizontal axis where the response time drops abruptly.  At those map scales the resampled pixel size is the same as the pixel sizes of the pyramids in the source data.  The processing time for resampling is the lowest at those scales because there is nearly a 1:1 match between source data pixels and resampled pixels.  Client applications which use this particular service will see dramatically faster response times if they are limited somehow to only those scales.  One way to do this is to use a tiled basemap layer.  Web mapping applications which use tiled basemaps are generally limited to only those map scales.  The most commonly used tiling scheme is the ArcGIS Online/Bing Maps/Google Maps tiling scheme (referred to hereafter as the “AGOL tiling scheme” for brevity).  The red-dotted vertical lines in the chart indicate the map scales for levels 7 – 12 of this tiling scheme.  Unfortunately those scales are not very close to the scales where this service performs it’s best.  There are two options for aligning source data pixels and tiling scheme scales:


    1. Build a custom basemap with a custom tiling scheme that matches the pixels sizes of the data.
    2. Sample or resample the data to a pixel size that matches the tiling scheme of the basemap.chart21.png

Chart 2. Performance of the binary overlay analysis service with different source data pixel sizes

The horizontal axis in chart 2 represents the "pixel size of the request" rather than map scale as in chart 1.  The orange graph shows the response times of another service configured identically to the first one in blue, except it uses source datasets that were up-sampled to 38 meter pixels using the Resample geoprocessing tool.  Up-sampling to 38 meters aligned the service’s fastest response times with the AGOL tiling scheme scales, which resulted in a significant decrease in processing time at those scales from approximately 1.5 seconds to about 0.5 seconds.  Furthermore, notice that performance is improved at nearly all scales except for the very largest.  This is most likely due to having all the source data at the same resolution (38m) instead of three (91.67m, 92.5m, 100m), and/or because the source data pixels are also aligned between datasets (accomplished by defining a common origin point for each resampled raster using the “Snap Raster” environment setting).


Admittedly, using the Resample tool to prepare data for analysis is not ideal because it results in second-generation data that is less accurate than the original.  This may be perfectly acceptable for applications intended to provide an initial survey-level analysis;however, it’s best to generate new first-generation data at the desired pixel size whenever possible.  For example, if you have access to land-class polygons, you could use them to generate a new first-generation raster dataset at the desired pixel size using the Polygon to Raster tool, rather than resampling an existing land-class raster dataset.


To determine by how much performance improved with 38 meter pixels, I calculated the percentage change in average response times for each scale and averaged the values over multiple scales.




Up-sampling the source data to 38 meter pixels reduced response times by  63.8% at the poorest-performing -target map scales!  38 meters was not my only option in this example.  I could have chosen a size that corresponded to one of the other tiling scheme scales.  The following table lists all the map scales of the AGOL tiling scheme, and the corresponding pixel sizes in units of meters, feet and Decimal Degrees.   The three columns on the right provide suggested values for sampling raster data.  These suggestions are not set in stone.  It’s not necessary to sample your data to exactly these recommended sizes.  The key is to choose a size that is slightly smaller than one of the sizes of the target tiling scheme scales.


Map Scales and Pixel Sizes for the ArcGIS Online/Bing Maps/Google Maps tiling scheme



By the way, matching the pixel sizes of your data with a basemap tiling scheme is also useful for workflows that involve static imagery overlaid onto a tiled basemap.  For those cases, you can build mosaic dataset overviews for viewing at smaller scales instead of raster pyramids.  One of the great things about mosaic dataset overviews is that you can define the base pixel size of overviews as well as the scale factor to match your target tiling scheme.  This way you don't have resample the source data to a new base pixel size in order to cater to any particular tiling scheme.

Resampling Method



The resampling method specified for an image service request also has an impact on performance.  The choice of which one to use should be based primarily on the type of data used in the analysis.  Chart 3 shows the performance of the binary overlay analysis service (with 38 meter data) with different resampling methods.


Chart 3. Response times of the binary overlay analysis service with different resampling methods




Bilinear resampling is the default method.  Here is how the response times for the other methods compared to bilinear averaged over the five map scales tested:


Raster Format


The storage format of the data can have a huge impact on performance.  For example, the response time of the binary overlay analysis service averaged over all map scales was 36% lower when the data was stored in the GeoTIFF format versus file geodatabase managed raster.  The Data Sources and Formats section of the  Image Management guide book recommends leaving the data in its original format unless it is in one of the slower-performing formats such as ASCII.  GeoTIFF with internal tiles is the recommended choice for reformatting because it provides fast access to the pixels for rectangular areas that cover only a subset of the entire file.Pixel Type and Compression


The pixel type determines the precision of the values stored in the data and can have a huge impact on performance.  In general, integer types are faster than floating-point types, and lower-precision types are faster than higher-precision types.  Compression of imagery can potentially increase or reduce performance depending on the situation.   For more information about the affect of compression on file size refer to the Image Management guide book section on Data Sources and Formats.  To assess the impact of pixel type and compression on the performance of data stored on a local hard drive, I tested a group of image services configured to perform an extremely intensive overlay analysis on 15 raster datasets.  The services were configured identically except for the pixel and compression types of the analysis data.  The tests were run at the map scale corresponding to the pixel size of the data.


Charts 5 & 6. Avg. response time and storage size vs. compression type for an image service that performs a complex overlay analysis



The following table shows the percentage change in response times with the reformatted datasets versus the original double-precision floating-point dataset.Table_pixeltype.png


On-the-fly Projection


On-the-fly projection is a very important feature of the ArcGIS platform.  It has saved GIS users like me countless hours of work by eliminating the need to ensure that every dataset in a map is stored in same coordinate system.  However, in some cases this flexibility and convenience may be costly when ultra-fast performance is required.   The following chart shows one of those cases.


Chart 8. Performance of a dynamic image service configured to perform a weighted overlay analysis.



Chart 8 shows the performance of a service which performs a weighted overlay analysis  on six datasets in an Albers projection.  The upper graph shows the performance when the output of the service is set to Web Mercator (Auxiliary Sphere).  The lower graph shows the performance when the output of the service is the same coordinate system as the data.  Performance without reprojection to Web Mercator improved by an average of 45% over all map scales.  This is a fairly extreme example.  The performance cost of reprojection is related to the mathematical complexity of the input and output projections.  Equal-area projections such as Albers are mathematically complex compared to cylindrical projections such as Mercator.  I have not run tests to prove this, but I expect that the performance cost of reprojection between two cylindrical projections such as UTM and Web Mercator would be less costly than seen in this example, and that a simple projection from geographic coordinates to Web Mercator would be even less costly.


To avoid on-the-fly projection you must ensure that all of your data is in the same coordinate system, including the basemap.  Most of the basemap services currently available from Esri on ArcGIS Online are in Web Mercator (auxiliary sphere).  So if you are going to use one of those basemaps, you would have to convert your data to the same coordinate system.  This can be an acceptable solution for some situations, but keep in mind that it results in second-generation data with less positional accuracy than the original source data.  Alternatively, you can create your own basemap in the same coordinate system as your data, and either publish it to an ArcGIS Server site or upload it to ArcGIS Online as a hosted map service.  If you take this approach, I recommend caching the basemap using a custom tiling scheme with scale levels that match the pixel sizes of your data.


Request Size


Request size is directly related to the size of the map window in the application and is specified in the REST API as the number of rows and columns of pixels in the output image.  To measure its impact on performance, I ran a series of tests at different request sizes on the weighted overlay analysis service that I used for the reprojection-on-the-fly tests.  I measured the average response times for request sizes ranging from 400x400 to 2200x2200, increasing at 100 pixel increments (e.g. 500x500, 600x600, etc…).  All of the tests were run at the map scale of 1:113386, which corresponds to the 30 meter pixel size of the source raster datasets.


Chart 9. Average response in MP/s for different request sizes for the weighted overlay service.


Chart 10. Average response time at different request sizes for the weighted overlay service.



Chart 9 shows that the throughput for this service levels off at a request size of approximately 1000x1000 pixels to about 1.5 – 1.6 MP/s.  Chart 10 shows that request size has a linear impact on performance.  This service is capable of providing sub-second response times for requests up to about 1,440,000 pixels, or a request size of 1200x1200.




Raster analysis can involve many stages of data processing and analysis.  Complex on-the-fly processing chains can place heavy processing loads on a server and contribute to sluggish performance.  Huge performance improvements can be achieved in some cases by pre-processing the data into a more efficient format for resampling and on-the-fly processing.


For applications which use tiled basemap layers, the greatest performance improvements are likely to be achieved by aligning the pixel sizes of the data with the scales of the basemap tiling scheme.  The section “Map Scales and Source Data Resolution” describes the theory behind this approach and provides a table with recommended pixel sizes for applications which use basemaps with the ArcGIS Online/Bing Maps/Google Maps tiling scheme.  Alternatively, developers can build basemaps with custom tiling schemes to align with the existing pixel sizes of the analysis data.


Another way to significantly reduce the processing load on a server in some cases is to avoid on-the-fly projection of the analysis data.  This is accomplished by ensuring that the basemap and the analysis data are in the same coordinate system.  The performance impact of on-the-fly projection varies depending on the input and output coordinate systems and is discussed in the section titled “On-the-fly Projection”.


The file format, pixel type, and compression type of the analysis data can also have a huge impact on performance.  GeoTIFF with internal tiles is recommended for situations where it’s necessary to re-format the data from a slower format.  Lower-precision pixel types give better performance than higher-precision types.  Pixel compression has the potential to either increase or decrease performance depending on the how the data is stored and accessed by the server.  These topics are discussed in the sections titled “Raster Format” and “Pixel Type and Compression”.


Client applications can also play a role in dynamic image service performance.  Service response times are the lowest when applications specify nearest neighbor resampling, followed by bilinear resampling.  And there is a direct relationship between service performance and the size of the map window in an application.  These topics are discussed in the sections titled “Resampling Method” and “Request Size”.

The Lab is pleased to announce the release of Raster Shader, a new addin for ArcMap 10.1!  Raster Shader is a tool for creating full-color, relief-shaded visualizations of single-band raster datasets.  Similar to the Image Analysis Window, Raster Shader leverages the power of raster functions to process the data dynamically as it is displayed in the map.  Since the processing occurs entirely in memory, it is extremely fast and does not affect the source datasets.


The idea for this tool came from projects we've done over the last couple of years involving large collections of single-band raster datasets containing climate data such as temperature averages and precipitation totals over time.  Our standard approach for visualizing these datasets is to use chained raster functions to apply a combination stretch/colormap renderer.  However, building and configuring raster function chains which produce an attractive and meaningful visual result can be a challenge if you are not familiar with the concepts and techniques involved.  This is where Raster Shader can help.  Not only does it build the raster function chain for you.  It also provides a simple interactive UI for adjusting key rendering parameters so you can see the results immediately in the map layer.  An interesting feature of Raster Shader is a histogram which shows how the colors in the colormap relate to the distribution of values in the data.  This relationship is based on stretch parameters which you can adjust with slider bars and/or text boxes.  Optionally, if the map contains a layer based on a DEM, you can add topographic context to the data with a hillshade effect and easily adjust parameters such as sun azimuth and sun angle.


Click here for the addin.

Filter Blog

By date: By tag: