Temporal data... statistical parameters over space
You have 12 months of data in some raster form... You want some statistical parameter... There are areas of nodata, the extents are all the same... and ... you have a gazillion of these to do.
Sounds like you have a 'cube'... the original 'space-time' cube' 
. You can pull space from a time slice... You can slice time through space. At every location on a 'grid' you have Z as a sequence over time.
Here is the code. I will use ascii files, but they don't have to be, you just need to prep your files before you use them.
Source data originally from .... here .... thanks Xander Bakker.
import os
import numpy as np
from textwrap import dedent, indent
import arcpy
arcpy.overwriteOutput = True
ft = {'bool': lambda x: repr(x.astype('int32')),
'float': '{: 0.3f}'.format}
np.set_printoptions(edgeitems=3, linewidth=80, precision=2, suppress=True,
threshold=80, formatter=ft)
np.ma.masked_print_option.set_display('-') # change to a single -
# ---- Processing temporal ascii files ----
# Header information
# ncols 720
# nrows 360
# xllcorner -180
# yllcorner -90
# cellsize 0.5
# NODATA_Value -9999
# ---- Begin the process ----
#
cols = 720
rows = 360
ll_corner = arcpy.Point(-180., -90.0) # to position the bottom left corner
dx_dy = 0.5
nodata = '-9999'
#
# ---- create the basic workspace parameters, create masked arrays ----
#
out_file = r'c:\Data\ascii_samples\avg_yr.tif'
folder = r'C:\Data\ascii_samples'
arcpy.env.workspace = folder
ascii_files = arcpy.ListFiles("*.asc")
a_s = [folder + '\{}'.format(i) for i in ascii_files]
arrays = []
for arr in a_s:
a = np.mafromtxt(arr, dtype='int', comments='#',
delimiter=' ', skip_header=6,
missing_values=nodata, usemask=True)
arrays.append(a)
#
# ---- A sample calculation from the inputs.... calculate the mean ----
#
N = len(arrays) # number of months... in this case
arr_shp = (N,) + arrays[0].shape # we want a (month, col, row) array
msk = arrays[0].mask # clone the mask... assume they are the same
e = np.zeros(arr_shp, 'int') # one way is create a zeros array and fill
for i in range(len(arrays)):
e[i] = arrays[i]
a = np.ma.array(e, mask=e*msk[np.newaxis, :, :]) # the empty array is ready
#
# ---- your need here... ie. Calculate a mean ----
#
avg = a.mean(axis=0) # calculate the average down through the months
#
# ---- send it out to be viewed in ArcMap or ArcGIS Pro ----
#
value_to_nodata = int(avg.get_fill_value())
out = avg.view(np.ndarray, fill_value=value_to_nodata)
g = arcpy.NumPyArrayToRaster(out, ll_corner, dx_dy, dx_dy)
g.save(out_file)
So the basic process is simple. I have coded this verbosely and used input parameters read manually from the ascii header since it is the quickest way.... and you would probably know what they are from the start.
So in this example... 12 months of some variable, averaged accounting for the nodata cells. Do the map stuff, define its projection... project it, give it some symbology and move on.
I will leave that for those that make maps.
Modify to suit... maybe I will get this into a toolbox someday
NOTE.....
Now in the linked example, there was a need to simply convert those to rasters from the input format. In that case you would simply consolidate the salient portions of the script as follows and create the output rasters within the masked array creation loop
......
for arr in a_s:
a = np.mafromtxt(arr, dtype='int', comments='#',
delimiter=' ', skip_header=6,
missing_values=nodata, usemask=True)
value_to_nodata = int(a.get_fill_value())
out = a.view(np.ndarray, fill_value=value_to_nodata)
r = arcpy.NumPyArrayToRaster(out, ll_corner, dx_dy, dx_dy)
out_file = arr.replace(".asc", ".tif")
r.save(out_file)
del r, a
.....So for either doing statistical calculations for temporal data, or for format conversion... there are options available where arcpy and numpy play nice.