Get annual maximum of daily rainfall and when it happen from time-series daily rainfall data in a year

readme.md

This is the revision from this version

The exisiting script worked perfectly on the small area of interest by reads all rasters into memory and then processes them. This can be very memory-intensive if we're dealing with large rasters (like global coverage).

Instead of loading all rasters into memory, we can process them one-by-one. This revision script now processes one raster at a time, updating the max value and its date, which should be more memory-efficient. We've added filename parsing to determine the date and then calculate the Julian day.

The script will also automatically detect the lower left point coordinates from the first raster in the directory and use them for the output rasters.

As not everyboody has ArcGIS license for arcpy, so I decided to create non-arcpy solution too.

max_and_date_arcpy.py

import sys
import numpy as np
import arcpy
import datetime

arcpy.env.overwriteOutput = True

# ---- You need to specify the folder where the tif files reside ----
src_flder = r"Y:\3days\2017"     # just change the year
out_flder = r"Y:\3days\2017_output"  # make a result folder to put stuff
out_year = src_flder.split("\\")[-1]
max_name = "{}\\idn_{}_max.tif".format(out_flder, out_year)
date_name = "{}\\idn_{}_day.tif".format(out_flder, out_year)

# Initialize arrays for max value and its corresponding date
first_raster_path = os.path.join(src_flder, arcpy.ListRasters()[0])
first_raster = arcpy.RasterToNumPyArray(first_raster_path, nodata_to_value=0.0)
max_values = np.zeros_like(first_raster, dtype=np.float32)
max_dates = np.zeros_like(first_raster, dtype=np.int32)

# Get the lower left point of the first raster
lower_left = arcpy.Point(arcpy.Raster(first_raster_path).extent.XMin, arcpy.Raster(first_raster_path).extent.YMin)

# Process each raster
arcpy.env.workspace = src_flder
rasters = arcpy.ListRasters()

for i in rasters:
    r = arcpy.RasterToNumPyArray(i, nodata_to_value=0.0)
    
    # Update max values and their corresponding dates
    condition = r > max_values
    max_values[condition] = r[condition]
    
    # Extract the date from the filename and convert to Julian day
    date_str = i.split('_')[-1].split('.')[0]
    date_obj = datetime.datetime.strptime(date_str, '%Y%m%d')
    julian_day = date_obj.timetuple().tm_yday
    max_dates[condition] = julian_day

# Save the results
out_max = arcpy.NumPyArrayToRaster(max_values, lower_left, x_cell_size=0.05, y_cell_size=0.05, value_to_nodata=0)
out_max.save(max_name)

out_when = arcpy.NumPyArrayToRaster(max_dates, lower_left, x_cell_size=0.05, y_cell_size=0.05, value_to_nodata=0)
out_when.save(date_name)

del max_values, max_dates, max_name, date_name, out_max, out_when

max_and_date_rasterio.py

import numpy as np
import rasterio
import os
import datetime

# ---- You need to specify the folder where the tif files reside ----
src_flder = r"Y:\3days\2017"     # just change the year
out_flder = r"Y:\3days\2017_output"  # make a result folder to put stuff
out_year = src_flder.split("\\")[-1]
max_name = "{}\\idn_{}_max.tif".format(out_flder, out_year)
date_name = "{}\\idn_{}_day.tif".format(out_flder, out_year)

# Initialize arrays for max value and its corresponding date
with rasterio.open(os.path.join(src_flder, os.listdir(src_flder)[0])) as src:
    meta = src.meta
    max_values = np.zeros_like(src.read(1), dtype=np.float32)
    max_dates = np.zeros_like(src.read(1), dtype=np.int32)

# Process each raster
for file_name in os.listdir(src_flder):
    if file_name.endswith('.tif'):
        with rasterio.open(os.path.join(src_flder, file_name)) as src:
            r = src.read(1)
            
            # Update max values and their corresponding dates
            condition = r > max_values
            max_values[condition] = r[condition]
            
            # Extract the date from the filename and convert to Julian day
            date_str = file_name.split('_')[-1].split('.')[0]
            date_obj = datetime.datetime.strptime(date_str, '%Y%m%d')
            julian_day = date_obj.timetuple().tm_yday
            max_dates[condition] = julian_day

# Save the results
with rasterio.open(max_name, 'w', **meta) as dst:
    dst.write(max_values, 1)

meta['dtype'] = 'int32'
with rasterio.open(date_name, 'w', **meta) as dst:
    dst.write(max_dates, 1)