I have 100 rasters of raw rainfall data in which the pixels each represent rainfall for that particular pixel on the landscape for a single year during the 100 year period. I have also created a raster that calculated the mean rainfall received in each pixel over those 100 years.

From these raw data rasters, I want to create a new set of 100 rasters (one for each year) that displays the standard deviation from the mean of each pixel (of the entire 100 years) for each pixel.

I have a feeling this will be likely solved with python coding, which I have a little bit of experience in, but I am unsure of what tool to use? It seems like Cell Statistics is on the right track, since it focuses on individual cells, but doesn't achieve what I am looking for.

Any other ideas on how to tackle this?

  • 1
    Because it is not possible that these rasters represent actual observations in every cell, they must have been interpolated from observations. As such, what you propose is going to reflect the interpolation method as much as--if not more than--actual patterns in the data and the resulting grid is going to be difficult if not impossible to interpret correctly. Why not instead analyze your monitoring stations separately? – whuber Nov 20 '13 at 16:32
  • You are correct that these data are interpolated from a series of datasets (I'm using PRISM rainfall data - prism.oregonstate.edu), so I don't have access to the individual monitoring stations. I also need continuous coverage for the analysis that I'm doing. I have checked the PRISM data against monitoring data that I do have, and they seem pretty good. – Colleen Nov 20 '13 at 17:52
  • 1
    Depending on the nature of the analysis, you could consider performing it on the stations and interpolating the results to get continuous coverage. You can recover the station data with reasonable accuracy by extracting the grid values at the station locations--so all you need is a dataset of their locations. – whuber Nov 20 '13 at 18:19

For each cell in your 100 rasters you know:

  • Current cell value (total rainfall per year?)
  • Mean rainfall per year at that cell

And you also know:

  • Standard deviation based on 100 years (population)

I would think the simplest approach would be to use the raster calculator in conjunction with a ModelBuilder iterator. You can write a short equation to place your current raster cell value about the mean and then assign the cell the appropriate value for standard deviation.

  • Thanks for the suggestion! I think I can play around and try this (although I'm not sure I know, but can calculate, the standard deviation based on 100 years). Basically, I want to create a raster, for each year, how much each cell deviates from the mean (of 100 years) for that same cell. – Colleen Nov 20 '13 at 17:53
  • 4
    Colleen, you must have used a local statistics calculation to obtain the mean at each cell. Another local stats calculation yields the SD at each cell. Your conversion is accomplished for each annual raster by subtracting that mean (which is a raster) and dividing the result by the SD (which is a raster, rather than a single number as intimated in this answer). That requires looping over the annual rasters but otherwise is simple and fast. – whuber Nov 20 '13 at 18:23
  • Yes, you are correct! Thanks for your patience with my trying to wrap my head around this. – Colleen Nov 20 '13 at 18:30
  • Agree with @whuber - working with a raster that represents SD makes more sense than dealing with a value in another (tabular/formulaic) format as I originally envisioned. – Radar Nov 20 '13 at 22:08

I assume that all your rasters are in the same folder, so you can lst them all with arccpy.ListRaster(). As you mentionned, you can then use cell statistics to get the mean, but also the standard deviation at each location

import arcpy
from arcpy import env
from arcpy.sa import *
env.workspace = "C:/workspace/with/data"
rasterList = arcpy.ListRasters("*","TIF")#listing all tifs in the workspace, could be another format
outSTD = CellStatistics(rasterList, "STD", "DATA") #note that "DATA" means  that the NoData values will not be used in the calculation.
outMEAN = CellStatistics(rasterList, "MEAN", "DATA")
for Individualraster in rasterList:
    #Compute how far you are from the mean value with respect to the STD. Under the assumption of a Gaussian distribution, 95% of your outputs will be between -1.96 and 1.96
    relativeSTD = (Raster(Individualraster) - outMEAN)/outSTD
    relativeSTD.save(Individualraster[:-4]+"relSTD.tif") #save under new name

Note that outMEAN and outSTD remain "inMemory" layers for "lazy" computation. If you want to store them, you need to save them in the script.

The ArcMap version is not specified: raster algebra could differ depending on the version.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.