Select to view content in your preferred language

Max value of time series raster and date information

4065
12
Jump to solution
07-29-2018 06:36 AM
BennyIstanto
Emerging Contributor

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‌

0 Kudos
12 Replies
DanPatterson_Retired
MVP Emeritus

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 ;)
BennyIstanto
Emerging Contributor

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 

Max in 2000

Date

0 Kudos
DanPatterson_Retired
MVP Emeritus

keep me posted if you need any more help Benny

0 Kudos