# How to calculate monthly maximum values from daily data in a 365 band grd?

I have a MINT2005.grd file which has 365 bands representing days of year from 2005. Each band is a 32x35 matrix representing the temperature at different latitude and longitude pair.

How do I calculate maximum value at each pixel position,so that the output contain 12 array of dimension 32x35 representing the month wise maximum value ?

• Is there sb. who knows : from osgeo import gdal ImportError: No module named osgeo why this error happened? – user49840 Mar 31 '15 at 14:44
• Hi sophia, could you please post this as a new question instead of posting it as an answer? You can always include the web address of this post in your question. – Joseph Mar 31 '15 at 15:05
• If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. – Dan C Mar 31 '15 at 15:20

You can use `np.max(axis=0)` to find the maximum along the z-axis of your array and get a 2d-array as a result. Then it is just a question of defining the monthly indices to loop through the original array.

``````import numpy as np
from osgeo import gdal

img = gdal.Open("/path/to/file")
dayspermonth = np.array([0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])

# instance empty array into which the results will be written
monthly_maxima = np.empty(shape=(12, 35, 32))

first = 0
last = 0
for i in range(len(dayspermonth)-1):  # loop through 12 months
# calculate first and last day of year for each month
first = first + dayspermonth[i]
last = last + dayspermonth[i+1]

# find maximum over z-axis, write into monthly layer
monthly_maxima[i, :, :] = np.max(input[first:last, :, :], axis=0)

# these are your monthly maxima at each grid-point
# Note: this does not account for leap years
january = monthly_maxima[0, :, :]
february = monthly_maxima[1, :, :]
# and so on...
``````

This can be easily done using GRASS GIS via "Processing" in QGIS 2: the r.series command of GRASS GIS accepts all 365 maps as input, then simply define method=maximum to filter the maximum for each pixel out of the stack. Adding also the "max_raster" method you can even determine which day of the year it was when reaching the respective maximum in a pixel stack (so: method=maximum,max_raster).