Automatic Raster Reclassification with ArcMap Model Builder Based on Raster Statistics

Alexandros Voukenas writes about the problem of automatic raster reclassification based on its statistics, using ArcMap’s Model Builder. This article offers a step-by-step guide with relevant screenshots and explanations, as well as a short overview of the theory behind this problem.

I’ve often considered starting a series of tutorials regarding common GIS and Remote Sensing problems that we face in various software and might be somewhat tricky to work around. Recently, I found what seemed to be a good start: Automatic raster reclassification with ArcMap’s Model builder based on the statistics of the raster.

In Remote Sensing, raster reclassification is a useful translation from continuous reflectance or index values to discrete classes that correspond to specific information. This is often performed in change detection studies, where the change of value between the same pixels in different dates corresponds to the level of change. For example, to study a burned area, the difference of Normalized Burn Ratio (NBR) is a good indicator for the burn severity:

ΔNBR= NBR(t1)-NBR(t2)

Obviously, a change value of 0 corresponds to an area that has not been burnt, and a value of 1 to an area that has been burnt severely. But how do we decide for all the in-between values? Traditionally, and considering that the pixel values in the image follow the normal distribution, the change severity is decided based on the raster statistics:

As we can see, in a population that follows (or approximates) the normal distribution, 64.2% of values range from (mean-1 standard deviation) to (mean+1 standard deviation). 27.2% of values range from (μ-2σ) to (μ-1σ) or from (μ+1σ) to (μ+2σ) etc. We can, therefore, translate this knowledge to classify the values of a change image into discrete values:

1. (μ-1σ to μ+1σ): No change
2. (μ-2σ to μ-1σ) or (μ+1σ to μ+2σ): Moderate change
3. (<μ-2σ) or (>μ+2σ): Severe change

This is all well-known standard Remote Sensing theory and the automatic reclassification of a raster is rather easy to pull off, using programming (with Python, R, or a similar language). Although I strongly encourage other young professionals like me, as well as students and young graduates, to start building confidence with programming as soon as possible and become independent from commercial software, I understand that not everyone is on the same level. Not to mention that often you will look for a quick an easy solution in a graphic environment and the time of building a custom script (which includes searching for documentation, debugging, etc.) does not seem “worth it”. In cases like these, graphic models of standard GIS software are good solution, and today we will see how to solve this problem easily with ArcMap’s Model Builder.

Suppose all we have as input data is one raster (for example, a change image) with continuous values and we want to automatically reclassify it into 3 classes based on its statistics: Mean value and standard deviation. To calculate these values, we will need the “Get Raster Properties” tool twice, once for mean and once for standard deviation:

Make sure to choose “MEAN” and “STD” in the “Property Type” tools. This will calculate mean and standard deviation, but here’s the tricky part (and mainly the reason I selected this topic today): Retrieving these values to use them in calculations for the reclassification. On GIS Stack exchangewe find a solution based on Raster Calculator, but this requires some specific syntax of variables and I was interested in a more simple and straightforward solution.

To get around this, we can use the “Create Constant Raster” Tool. As the title suggest, it will create a raster with a constant value and we can set these values to be the mean and standard deviation. We also have to select the output data type as parameter for Create Constant Raster:

To complete the tool, we have to set spatial reference as parameters: Cell size and Output extend. We can do that by connecting the output raster to the tool directly twice, once for each parameter:

In addition to that, we can keep only one “Output data type” parameter and connect it to both Create Constant Raster:

So, with this indirect method, we have retrieved the statistical values and we’re in position to make calculations with them. We have to make the calculations:

1. (μ-1σ to μ+1σ)
2. (μ-2σ toμ-1σ) OR(μ+1σ toμ+2σ)
3. (<μ-2σ) OR (>μ+2σ)

We first use the Times tool for standard deviation to multiply it by 2. Then, with the Plus and Minus tools and the outputs “STDV Raster” and “Mean Raster” we can perform all the calculations. Let’s start with the first interval:

Excuse me if some operators get mixed up, this is because I constantly use the “Auto layout” option of Model Builder. So, we have created the rasters μ+1σ and μ-1σ. And we want all the values of our initial raster within that interval, to be part of the same class. We can easily achieve that with a combination of order and Boolean tools, such as the “Less than”, “Greater than”, “Boolean OR” and “Boolean AND”:

So, what we did was quite simple. We calculated two rasters: (μ+1σ) and (μ-1σ). If the values of our initial raster ranged within these two values (“Less than Equal” and “Greater than Equal” tools), then a binary raster is created (TRUE if they fall within, FALSE otherwise-”Boolean AND” operator).

We can create the other two classes similarly using the 2x STDV Raster and the same tools, so let me fast forward to the completed 2nd class:

I’ve purposely cut out the rest of the model which calculates class 1 in order to not make it look too complicated. The second class consists of two separate subclasses (what I call “Upper” and “Lower”, or you can call them “Left” and “Right”) so we need an additional “Boolean OR” Tool to make the final Class 2: Moderate change.

The final class is created similarly, based only on the inequality tools with μ+2σ and μ-2σ:

Now, we have 3 separate TRUE/FALSE rasters. If you look closely to the inequality tools, you will notice that I’ve picked them carefully in order to ensure non-overlapping intervals:

1. Class 1: [μ-1σ, μ+1σ]
2. Class2: [μ-2σ, μ-1σ) or(μ+1σ, μ+2σ]
3. Class 3: [min, μ-2σ) or (μ+2σ, max]

So what we need for our final thematic layer is to “merge” together these 3 rasters. To do that, we can easily convert the binary TRUE/FALSE to Integer 1/0. Then, we multiply each Integer raster to another integer in order to tell them apart and finally add them all together, and each discrete value (0, 1 or 2) will correspond to that specific class:

I mentioned earlier that we might prefer a nice graphic solution instead of writing code. Truth is, when you look at the whole model, it doesn’t look so elegant:

But the truth is, the steps are taken are quite easy and the flow is understandable. Also, when I worked with this problem, I was working only on an initial binary classification, for example:

<μ+2σ: 1, ELSE 0

This needs a fairly smaller model, but it’s just as easy. Nevertheless, we have created a ton of intermediate rasters, so make sure they ARE intermediate, by right clicking on them:

This way, using Model Delete Intermediate Variables, they will be deleted from our database. One thing left to do: Run the model and see the result:

We have our three classes, as expected. The processing time was around 4 minutes in an 8GB, i5 machine, which is not bad considering the original file came from an ENTIRE Sentinel 2 scene (290 x 290 km, ~500 MB).

This was the first GIS problem that I decided to present and I’m hoping to establish a series of similar problems. I hope it was a useful solution and I’m really looking forward to feedback regarding clarifications, possible improvements, or proposals for other problems to tackle!

Bibliography

ESRI, n.d. What is Model Builder?. [Online]
Available at: http://desktop.arcgis.com/en/arcmap/10.3/analyze/modelbuilder/what-is-modelbuilder.htm
[Accessed February 2019].

Lunetta, R. S., 1999. Remote sensing change detection: environmental monitoring methods and applications. London: Taylor and Francis.

About the Author

Alexandros Voukenas received his diploma (MSc equivalent) in Surveying Engineering from Aristotle University of Thessaloniki in 2017 and has been working as a GIS & Remote Sensing Engineer in the private sector ever since. His main areas of occupation have been the production of reference mapping using Geographic Information Systems and satellite imagery, and to perform feasibility research studies on remote sensing projects.

Voukenas is passionate about solving of real life problems and he has won a prize in four innovation hackathons, including the Fabspace 2.0 Greece (1st place), the NASA Space Apps Challenge Thessaloniki (2nd place) and the Crowdhackathon smarticity (4th place). He is also an active researcher on the topics of urban data analysis and remote sensing. His other research interests are spatial analysis & statistics, GIS, geospatial technologies, photogrammetry and computer vision.

Alexandros Voukenas can be contacted via email.

More from Alexandros Voukenas:

• Calculating Land Use Mix with GIS

Related

• How to Use ArcGIS Pro to Map Flood Susceptibility
• How to Use ArcGIS Pro to Map Urban Heat Islands
• How to Use ArcGIS Pro to Estimate Areas Prone to Malaria