Accepted Solutions
Benny... do a manual check.
Set the color ramp for the max to 'Precipitation'
The other raster is the 'day' that the maximum occurred.
Let me know... it takes less than 30 seconds to do a run... 30 * 30 is less time to automate it to do
everything 
you just need make folders in line 44, 45... put all the raster folders into line 44's folder and the results have to go into a separate folder
# -*- coding: utf-8 -*-
"""
climate
=======
Script : climate.py
Author : Dan_Patterson@carleton.ca
Modified : 2018-07-30
Purpose:
Process daily climate data... for this time, the daily max by
location. Bali, is the locale.
Useage :
References
----------
`<https://community.esri.com/thread/218763-max-value-of-time-series-raster-
and-date-information#comment-788295>`_.
`<http://www.bali-paradise.com/news/weather/`>_.
`<http://pro.arcgis.com/en/pro-app/arcpy/data-access/numpyarraytotable.htm>`_.
`<http://pro.arcgis.com/en/pro-app/arcpy/data-access/tabletonumpyarray.htm>`_.
---------------------------------------------------------------------
"""
import sys
import numpy as np
import arcpy
arcpy.env.overwriteOutput = True
ft = {'bool': lambda x: repr(x.astype(np.int32)),
'float_kind': '{: 0.3f}'.format}
np.set_printoptions(edgeitems=5, linewidth=120, precision=2, suppress=True,
threshold=100, formatter=ft)
np.ma.masked_print_option.set_display('-') # change to a single -
script = sys.argv[0] # print this should you need to locate the script
# ---- You need to specify the folder where the tif files reside ------------
#
src_flder = r"C:\Temp\climate\1981" # just change the year
out_flder = r"C:\Temp\climate\Results" # make a result folder to put stuff
out_year = src_flder.split("\\")[-1]
max_name = "{}\\Bali_{}max.tif".format(out_flder, out_year)
date_name = "{}\\Bali_{}_day.tif".format(out_flder, out_year)
# ---- Process time ----
arcpy.env.workspace = src_flder
rasters = arcpy.ListRasters()
arrs = []
for i in rasters:
r = arcpy.RasterToNumPyArray(i, nodata_to_value=0.0)
arrs.append(r)
# ---- Convert to an integer dtype since the float precision is a bit much
a = np.array(arrs)
m = np.where(a > 0., 1, 0)
a_m = a * m
a_m = np.ndarray.astype(a_m, np.int32)
a_max = np.max(a_m, axis=0)
a_w = np.argmax(a_m, axis=0)
a_when = np.ndarray.astype(a_w, np.int32)
out_max = arcpy.NumPyArrayToRaster(a_max,
arcpy.Point(114.4317366,-8.8492618),
x_cell_size=0.049242234706753,
y_cell_size=0.049242234706753,
value_to_nodata=0)
out_max.save(max_name)
out_when = arcpy.NumPyArrayToRaster(a_when,
arcpy.Point(114.4317366,-8.8492618),
x_cell_size=0.049242234706753,
y_cell_size=0.049242234706753,
value_to_nodata=0)
out_when.save(date_name)
del a, m, a_m, a_max, a_w, a_when, max_name, date_name, out_max, out_when
# ---- repeat ad nauseum since this is free ;)Hi Dan, thank you for reply
I am not familiar with numpy and friend, I can only use script, without knowledge write a code. You and Xander Bakker always helping me solve my previous problem 
My latest question to you and Xander about "percentile from raster timeseries".
From your suggestion, which tool should I use to compare? Is it doable using Model Builder raster iteration?
Benny
ahhh yes How to calculate the percentile for each cell from timeseries raster
# ---- lines 19 and 20
a_m = np.array(arrs_m)
a_max = np.nanmax(a_m, axis=0)
# ------ now an example --------------------
a = np.random.randint(1, 10, size=(4, 5, 6))
a # ---- years
array([[[6, 4, 7, 1, 1, 1],
[1, 1, 3, 9, 6, 7],
[1, 6, 9, 9, 3, 6],
[6, 8, 3, 5, 8, 4],
[5, 8, 7, 6, 1, 5]],
[[9, 1, 6, 7, 9, 1],
[1, 3, 7, 8, 2, 9],
[2, 1, 8, 7, 5, 3],
[7, 9, 1, 3, 6, 4],
[1, 2, 1, 2, 2, 1]],
[[9, 9, 3, 4, 3, 6],
[2, 3, 9, 6, 5, 1],
[7, 8, 8, 5, 1, 7],
[6, 7, 2, 8, 4, 6],
[1, 8, 1, 8, 6, 2]],
[[3, 4, 5, 2, 6, 1],
[6, 3, 2, 6, 9, 2],
[3, 4, 9, 3, 7, 7],
[6, 4, 4, 9, 6, 3],
[6, 9, 2, 2, 4, 8]]])
a_max = np.nanmax(a, axis=0)
a_max # ---- The maximum by location regardless of year
array([[9, 9, 7, 7, 9, 6],
[6, 3, 9, 9, 9, 9],
[7, 8, 9, 9, 7, 7],
[7, 9, 4, 9, 8, 6],
[6, 9, 7, 8, 6, 8]])
# ---- now the 'time' that the max occurred
np.argmax(a, axis=0)
Out[8]:
array([[1, 2, 0, 1, 1, 2],
[3, 1, 2, 0, 3, 1],
[2, 2, 0, 0, 3, 2],
[1, 1, 3, 3, 0, 2],
[3, 3, 0, 2, 2, 3]], dtype=int64)Now this requires interpretation...
let's look at line 37 row 0
9, 9, 7, 7, 9, 6
now line 47 row 0
1, 2, 0, 1, 1, 2
the max occurred at those locations at.
For example top left corner, max is 9 and the time period was 1 (counting from 0 )
next over, max is 9, time period was 2
next over, max is 7, time first period 0
next over, max is 7, time period 1
next over, max is 9, time period 1
and the last, max is 6, time period 2
repeat.
Now if you want to find out which time period had the most maximums by location... I can do that.
If you want to see if there is any clustering spatially in the time of maximum... ditto
how about …. fill in the blank …..
So see if you can resurrect the thread and give it whirl yourself
Benny … 1981 max …
Looks ok?
I will spruce itup after you give the ok.
I had to make some minor changes to the code to simplify a few things.
Now... do you really need the precision to 1/100 of a millimeter? 
or would the nearest mm be ok?
