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

DanPatterson_Retired

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

(1)  Reading the raster files

<SPAN class="keyword token">import</SPAN> numpy <SPAN class="keyword token">as</SPAN> np
<SPAN class="keyword token">import</SPAN> arcpy

arrs <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">]</SPAN>
folder <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"C:\Data\rasters"</SPAN>&nbsp; <SPAN class="comment token"># ---- specify a folder containing all the rasters ----</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>workspace <SPAN class="operator token">=</SPAN> folder
rasters <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListRasters<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"*"</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"TIF"</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">for</SPAN> raster <SPAN class="keyword token">in</SPAN> rasters<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; arrs<SPAN class="punctuation token">.</SPAN>append<SPAN class="punctuation token">(</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>RasterToNumPyArray<SPAN class="punctuation token">(</SPAN>raster<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>

a <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>array<SPAN class="punctuation token">(</SPAN>arrs<SPAN class="punctuation token">)</SPAN> <SPAN class="comment token"># ---- the master array ----‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

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 ----‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

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 <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>Raster<SPAN class="punctuation token">(</SPAN>r01<SPAN class="punctuation token">)</SPAN>&nbsp; <SPAN class="comment token"># ---- read the first raster we loaded ----</SPAN>

rast <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>Raster<SPAN class="punctuation token">(</SPAN>r01<SPAN class="punctuation token">)</SPAN>&nbsp; <SPAN class="comment token"># ---- simple </SPAN>

dir<SPAN class="punctuation token">(</SPAN>rast<SPAN class="punctuation token">)</SPAN>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># ---- get raster information... some snipping here ----</SPAN>

<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">.</SPAN><SPAN class="punctuation token">.</SPAN><SPAN class="punctuation token">.</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'bandCount'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'catalogPath'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'compressionType'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'extent'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'format'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'hasRAT'</SPAN><SPAN class="punctuation token">,</SPAN>
 <SPAN class="string token">'height'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'isInteger'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'isTemporary'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'maximum'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'mean'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'meanCellHeight'</SPAN><SPAN class="punctuation token">,</SPAN>
 <SPAN class="string token">'meanCellWidth'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'minimum'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'name'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'noDataValue'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'path'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'pixelType'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'save'</SPAN><SPAN class="punctuation token">,</SPAN>
 <SPAN class="string token">'spatialReference'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'standardDeviation'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'uncompressedSize'</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'width'</SPAN><SPAN class="punctuation token">]</SPAN>

<SPAN class="comment token"># ---- Save the result out to a new raster ------</SPAN>

r01 <SPAN class="operator token">=</SPAN> rasters<SPAN class="punctuation token">[</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">]</SPAN>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># --- rasters have the same extent and cell size</SPAN>
rast <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>Raster<SPAN class="punctuation token">(</SPAN>r01<SPAN class="punctuation token">)</SPAN>

lower_left <SPAN class="operator token">=</SPAN> rast<SPAN class="punctuation token">.</SPAN>extent<SPAN class="punctuation token">.</SPAN>lowerLeft&nbsp; <SPAN class="comment token"># ---- needed to produce output</SPAN>
cell_size <SPAN class="operator token">=</SPAN> rast<SPAN class="punctuation token">.</SPAN>meanCellHeight&nbsp;&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># ---- we will use this for x and y</SPAN>

f <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"c:\temp\result.tif"</SPAN>
r <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>NumPyArrayToRaster<SPAN class="punctuation token">(</SPAN>a<SPAN class="punctuation token">,</SPAN> lower_left<SPAN class="punctuation token">,</SPAN> cell_size<SPAN class="punctuation token">,</SPAN> cell_size<SPAN class="punctuation token">)</SPAN>
r<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN>f<SPAN class="punctuation token">)</SPAN>
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

(4) Saving and reloading the arrays

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<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN>r<SPAN class="string token">"c:\temp\myarray.npy"</SPAN><SPAN class="punctuation token">,</SPAN> a<SPAN class="punctuation token">)</SPAN>&nbsp; <SPAN class="comment token"># ---- 'a' is the multi dimensional array</SPAN>

a<SPAN class="punctuation token">.</SPAN>shape&nbsp; <SPAN class="comment token"># ---- 31 days, 100 x 100 cells </SPAN>
<SPAN class="punctuation token">(</SPAN><SPAN class="number token">31</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">100</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">100</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token"># ---- reload it </SPAN>

im_back <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>load<SPAN class="punctuation token">(</SPAN>r<SPAN class="string token">"c:\temp\myarray.npy"</SPAN><SPAN class="punctuation token">)</SPAN>

im_back<SPAN class="punctuation token">.</SPAN>shape
<SPAN class="punctuation token">(</SPAN><SPAN class="number token">31</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">100</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">100</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token"># ---- check to see if it is good</SPAN>

np<SPAN class="punctuation token">.</SPAN>all<SPAN class="punctuation token">(</SPAN>a <SPAN class="operator token">==</SPAN> im_back<SPAN class="punctuation token">)</SPAN>&nbsp; <SPAN class="comment token"># ---- some numpy magic ;)</SPAN>
<SPAN class="token boolean">True</SPAN>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

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.

Array_small_31.zip

DanPatterson_Retired

Luke... you can tile your data prior to running or there is a 'tile' option (2nd code example in RasterToNumPy array )) that one could employ but I am usually on using rasters < 5k x 5k ish cell wise. 

Most statistics can determined by either tiling in the x-y plane or by tiling the z (time) plane... reassembling the result at the edge (ie.. batch a year's data by month rather than doing the whole year at once)

Nan in most circumstances is a great way of just using the nan-statistics option and the choice between nan or a masked array often depends on the statistic being calculated.  The best example is the fairly recent addition of nan.percentile to the compliment.

'Slow' is a matter of perception sometimes so I haven't done much testing with that end of things compared to the actual Cell statistics tool.  Other options depend on whether you want to calculate the statistic directly (ie nan.std) or do it in steps.  If 0 (zero) is not present as a valid observation, then it can be used to generate a final mask.

Where it gets tricky is when the nodata values aren't in the same location over time (ie truly missed observations) then your only options are nan statistics or the masked array.