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")
input = img.ReadAsArray()
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).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.