0

I have been trying to read in a stack of images for a timeseries and created an output raster that contains the max pixel values of the timeseries.

The arrays were created using numpy

img_stack = np.empty((image_width, image_length, len(image_files)), np.dtype('f'))      
final_values = np.empty((image_length, image_width), np.dtype('f'))

I imported the image using the following loop:

for i, fname in enumerate(image_files):
    img = gdal.Open(os.path.join(fn + ("\{0}".format(fname)))).ReadAsArray()
    img_stack[:, :, i] = img

then I read the values into a 2D list using the following nested loop:

for x in range(image_length):
    for y in range(image_width):
        final_values[x][y] = 0
        for z in range(len(image_files)):
            if (final_values[x][y] < img_stack[y, 0, z]):
                final_values[x][y] = (img_stack[y, 0, z])
            else:
                next

then use

output_raster = dr.Create(final, image_width, image_length, 1, gdal.GDT_Float32)
output_raster.GetRasterBand(1).WriteArray(final_values)

the output contains the correct pixel values, but is rotated 90 degrees to the left.

I believe the issue is due to GDAL reading the image as an array first by width and then by length, but I am not sure.

Is there an issues with my code or is there a function I could use to rotate it to the proper orientation?

Thanks

  • 1
    Is the output truly rotated or has it actually been reflected? A reflection would result from transposing the row and column coordinates whereas a rotation would suggest the row and column coordinates have been transposed and the order of the rows has been reversed. (These are the two details that cause the most problems when converting among various raster and array indexing conventions.) – whuber Apr 3 '14 at 22:51
  • 1
    the output is neither reflected nor rotated, it's gibberish. As noted in Marc Pfister's reply, the code is checking img_stack[y, 0, z], but setting final_values[x,y]. The x var indexes image_length, which makes the code somewhat inscrutable to the casual reader, but also means that for a non-square array, img_stack should be getting array bound errors if width > length. In addition, only the first (0) column values are read, but written into the max for every column (though transposed in the process) – Llaves Apr 4 '14 at 5:10
  • 2
    another thing - this approach will eat memory like crazy if the images are large and there are many images. It would be more space efficient to read one image and compare its pixel values to the current high value, then read the next image, and so on. – Llaves Apr 4 '14 at 5:12
1

I'm not sure why you're accessing

img_stack[y, 0, z]

instead of

img_stack[x, y, z]

But let's make it easier - if you want the maximum values in a numpy array along your z axis, you should be able to use:

final_values = numpy.amax(img_stack, axis=2)

see numpy.amax

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.