Hi,
I have daily rainfall data in a year in geotiff format with file naming convention "rainfall_yyyymmdd.tif" and I would like to get maximum rainfall value per year per pixel and information about which date the max value came from.
Getting the max value is doable using Cell Statistics, but how to get information about date of the max value?
I am planning to do the calculation for 36 years data (1981 - 2017).
Thank you, Benny
raster timeseries arcpy python cellstatistics precipitation rainfall
Solved! Go to Solution.
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 ;)
Awesome!
Yes, I confirmed it takes less than 30 seconds.
I tested using Indonesia data (whole country), and modified line 67 and 74. I guess its the lower left coordinates.
Thank you very much Dan Patterson you made my day.
FYI, I will use this script to see the relation between max rain+date, flood frequency and extreme rainfall more than threshold (generated using recurrence interval).
This question is still related with my previous question about raster percentile from timeseries data
keep me posted if you need any more help Benny