Raster change calculations and zonal statistics at scale
Raster data is often a snapshot in time, for instance representing land cover, vegetation, or impervious surfaces. But the real value shows up when you compare those snapshots.
With the raster capabilities in ArcGIS GeoAnalytics Engine, you can compute and summarize change at scale directly in your Spark workflows, no need to step outside your pipeline in your analytics platform. In this post, we’ll walk through a simple but powerful pattern: calculating land use / land cover (LULC) change and summarizing it with zonal statistics.
This type of workflow has incredible value across domains, for instance:
- In retail and commercial planning, you can track how development is expanding around store locations—helping identify where demand is growing before it shows up in sales data.
- In insurance, land cover change can act as an early signal for risk. Urban expansion into wildfire-prone areas or vegetation loss in flood zones can materially impact underwriting decisions.
- In energy and natural resources, it’s a straightforward way to monitor vegetation loss near infrastructure, track reclamation progress, or quantify environmental impact within regulated areas.
- In urban planning, it’s one of the simplest ways to measure sprawl, loss of green space, or changes in land use policy outcomes.
Let’s look at an example of how we can do these types of analysis with GeoAnalytics Engine working in an Apache Spark distributed computing environment.
Calculating change
We’ll begin with two land use / land cover (LULC) raster datasets from different time periods. In this example we will use LULC data from the Multi-Resolution Land Characteristics Consortium. Specifically, we will use Land Cover data from 2020 and 2024. These are US national-scale datasets with a 30-meter resolution.
Each of the datasets provides a land cover value for every 30-meter cell across the conterminous United States. As an example, let’s take a look at the 2020 data:
The land use classes are divided into the following categories:
To calculate change with GeoAnalytics Engine, we can make a composite raster with the two years of data using the CompositeBands tool, and then use RT_Apply to apply a change calculation. In this case, we are just looking for all cells where the land use value differs between the two years. If they differ, we assign the value from the most recent year, otherwise we return no value.
from geoanalytics.tools import CompositeBands
# make compsite with both years of landcover data in separate bands
composite = CompositeBands()\
.addRaster(landcover_2020, "raster", [1])\
.addRaster(landcover_2024, "raster", [1])\
.run()
change_lambda = lambda px: F.when(px.values[1] != px.values[0], px.values[1]).otherwise(None)
luc_change = composite.select(RT.select_bands(RT.apply("composite_raster", 2, change_lambda), [2]).alias("change"))This gives us a new calculated raster with the land use change. We can zoom in to see how this looks for an area near San Diego, California – the individual raster cells with changes in land use between the two years are shown:
We could also process the data further to look at the change as a matrix to see how the land use classes have changed across years:
Summarizing change where it matters
Once you have a change raster, the next step is usually: where does this change matter?
That’s where zonal statistics comes in, summarizing raster values inside polygons like counties, trade areas, or risk zones.
This gives you a quick way to quantify how much change occurred, and how complex it was, within each zone. As an example, we will use hexagonal bins to analyze change across our area of interest near San Diego.
from geoanalytics.tools import ZonalStatistics
result = ZonalStatistics() \
.setZones(bin_geometries) \
.setZoneIdColumns("bins") \
.setRasterColumn("change") \
.includeZoneGeometry(True)\
.run(luc_change)Now, for each of our hexagonal bins we know the details on all of the cells changing value in that zone – the count, majority and minority class, percent covered by that class, etc.
Focusing on the changes you care about
Often, you don’t need all transitions, you care about specific ones. For example, exploring all of the changes into a “Developed” land use class. Note that these could be changes from undeveloped to developed or between intensities of development. We can explore these patterns with a filter to focus the view on just the regions where the majority change is to become a Developed land use class and where the change percent is greater than 50%.
Now we can more easily see the regions where the development is focused, and where the intensity of development is increasing the most.
A simple pattern that scales
Change detection is a foundational raster workflow, but what’s new here is how easily it scales. While our examples in this post demonstrated calculations for a small part of the United States near San Diego, California, you can run this across large areas, integrate it with feature data, and push results straight into downstream analytics.
Scaling this with Spark to run across three 16-core machines, we can calculate land use change using the 30 meter national scale raster datasets and then compute the zonal statistics for ~20 million level 8 H3 bins (~0.78 square kilometers) in about 7 minutes!
If you’re already using GeoAnalytics Engine, this is a natural extension of what you’re doing today, now including raster data in the mix.
Conclusion
Raster support in GeoAnalytics Engine opens up a new class of workflows, bringing large-scale surface analysis directly into your distributed pipelines alongside vector data. Whether you’re detecting land cover change, monitoring environmental conditions, or tracking development patterns, these capabilities make it much easier to move from raw raster data to actionable insight.
We’re excited to see how you apply these patterns in your own workflows, especially as you combine raster and feature analytics in new ways. Let us know what you’re building and where raster is helping unlock new value!
