# Cell Statistics made easy... raster data over time... the "Special Analyst"

Blog Post created by Dan_Patterson on Feb 6, 2018

This comes up time and again. ....I have 'X' rasters for an area and I need to determine Y over time....

Here is an example.... ....calculating percentiles....

Cell Statistics... is the obvious choice.  But you don't have the Spatial Analyst extension? Doesn't matter.  All the tools that you need are already provided if you are using current versions of ArcMap or ArcGIS PRO.

Do some work

``import numpy as npimport arcpyarrs = []folder = r"C:\Data\rasters"  # ---- specify a folder containing all the rasters ----arcpy.env.workspace = folderrasters = arcpy.ListRasters("*", "TIF")for raster in rasters:    arrs.append(arcpy.RasterToNumPyArray(raster))a = np.array(arrs) # ---- the master array ----``

Once the rasters are all converted to arrays, they were appended to a list.  A master array was created in line 9.

As an example, 31 rasters in tif format were read in.  These are small rasters, but they can be as large as you can process.  You are now ready to calculate values for them using numpy.

(2)  Getting the statistics

``median = np.median(a, axis=0) # ---- pick a statistic ---- array([[ 59.,  54.,  36., ...,  57.,  46.,  46.],       [ 43.,  45.,  59., ...,  38.,  51.,  51.],       [ 35.,  50.,  57., ...,  47.,  55.,  65.],       ...,        [ 43.,  52.,  40., ...,  56.,  62.,  60.],       [ 45.,  57.,  39., ...,  45.,  44.,  48.],       [ 49.,  57.,  50., ...,  56.,  50.,  58.]]mins = np.min(a, axis=0)array([[1, 3, 0, ..., 6, 4, 5],       [1, 5, 6, ..., 6, 7, 2],       [3, 2, 4, ..., 2, 5, 4],       ...,        [4, 0, 1, ..., 3, 4, 6],       [0, 0, 0, ..., 5, 1, 0],       [4, 1, 1, ..., 0, 3, 7]])# ---- name a statistic, with or without nodata values, it can be done ----``

Now when you are done your calculations you will probably want to make a raster so you can bask in the beauty of your efforts.

(3)  Reading the raster information and saving the result

``rast = arcpy.Raster(r01)  # ---- read the first raster we loaded ----rast = arcpy.Raster(r01)  # ---- simple dir(rast)                 # ---- get raster information... some snipping here ----[..., 'bandCount', 'catalogPath', 'compressionType', 'extent', 'format', 'hasRAT', 'height', 'isInteger', 'isTemporary', 'maximum', 'mean', 'meanCellHeight', 'meanCellWidth', 'minimum', 'name', 'noDataValue', 'path', 'pixelType', 'save', 'spatialReference', 'standardDeviation', 'uncompressedSize', 'width']# ---- Save the result out to a new raster ------r01 = rasters[1]                    # --- rasters have the same extent and cell sizerast = arcpy.Raster(r01)lower_left = rast.extent.lowerLeft  # ---- needed to produce outputcell_size = rast.meanCellHeight     # ---- we will use this for x and yf = r"c:\temp\result.tif"r = arcpy.NumPyArrayToRaster(a, lower_left, cell_size, cell_size)r.save(f)``

Now, the nice thing about arrays is that you can save a 3D to a file for reloading later (or any dimension for that matter).

If you have a 3D array like we have, this is kind of like a big 'zip'.

``np.save(r"c:\temp\myarray.npy", a)  # ---- 'a' is the multi dimensional arraya.shape  # ---- 31 days, 100 x 100 cells (31, 100, 100)# ---- reload it im_back = np.load(r"c:\temp\myarray.npy")im_back.shape(31, 100, 100)# ---- check to see if it is goodnp.all(a == im_back)  # ---- some numpy magic ;)True``

That's about all.  Not covered here is nodata cells (easy, you can assign one) and statistics for nodata (easy, every stat has a nodata equivalent).

Try your skills with the attached *.npy file.