Averaging image/bands on a per-pixel basis without having to explicitly name bands?

[Thiago Sanna Freire Silva]

Hello everyone,

I have several single-band raster files which I would like to average on a per-pixel basis. I though of using imagecalc or gdal_calc, but they both assume you name each band and write the expression explicitly. I'd rather like to to something like this:

import glob

import rsgislib

imgList = glob.glob("*_HH.tif")

outImg = "average_HH.tif"

# imaginary function 

applyfunctionperpixel(imgList,outImg, function=mean, "GTiff", rsgislib.TYPE_32UINT)

Any suggestions?

Thanks,

Thiago

[Nathan Thomas]

Hi Thiago,

I do not think there is currently a function in RSGISLib that will do this for you. Below is some python that will achieve what you want, although it may be a little ugly.

import rsgislib

from rsgislib import imagecalc

from rsgislib.imagecalc import BandDefn

import os

outputImage = ‘outputfile.kea'

inDir = ‘input_path'

# Gather files

files = [file for file in os.listdir(inDir) if file.endswith('.kea')]

bandDefns = []

explist = []

for file in files:

    # Define each single file band

    bandDefns.append(BandDefn('b' + str(files.index(file)), os.path.join(inDir, file), 1))

    # Make a list of the file number (i.e. [1,2,3,4,5...])

    explist.append('b' + str(files.index(file)))

gdalformat = 'KEA'

datatype = rsgislib.TYPE_32FLOAT

# Create the expression by joining the band bumber to a 'b' for each file

# and dividing by the total number of files

expression = str(' + ').join(explist).join('()') + '/' + str(len(explist))

imagecalc.bandMath(outputImage, expression, gdalformat, datatype, bandDefns)

This worked for a quick example that I made for myself, but let me know if you get any problems. You will still have to write out the expression by hand, which may be difficult if you want to do more complex calculations other than sum, average etc, but this should do what you need for now.

Thanks, Nathan

Nathan M Thomas

CalTech Postdoc at JPL

Radar Science and Engineering

4800 Oak Grove Dr.

300-356A

Pasadena, CA 91109

Phone: +1 818 354 4681

--
You received this message because you are subscribed to the Google Groups "RSGISLib Support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-suppo...@googlegroups.com.
To post to this group, send email to rsgislib...@googlegroups.com.
Visit this group at https://groups.google.com/group/rsgislib-support.
For more options, visit https://groups.google.com/d/optout.

[Daniel Clewley]

Hi Thiago,

There is the stackStats (http://rsgislib.org/rsgislib_imageutils.html#rsgislib.imageutils.stackStats) function in RSGISLib which I think does what you need, but all the files need to stacked into a single file, you can use something like:

gdalbuildvrt -separate out_stack.vrt In_files*.tif

to produce a virtual stack and then use this as an input. I haven’t used this function in a while so let me know if there are any problems with it.

Another option is to use RIOS, I tend to use this for applying more complicated functions to images. It works using numpy arrays so you can use a lot of the built in numpy functions

I posted an example on stackexchange a while ago in response to a similar question:

https://gis.stackexchange.com/questions/88996/per-pixel-statistical-calculations-on-a-raster-stack-using-gdal/107344#107344

You should be able to adjust this code do what you need.

Thanks,

Dan

[Nathan Thomas]

Hi Dan,

I originally thought of using RIOS too, but I wasn’t sure how to handle multiple files. Is there a way to pass a directory/multiple files into RIOS? Ive only used it to pass in single str or int arguments. Thus, I went with the hacked together bandmath approach.

Otherwise, stackStats is the best way to go.

Cheers, Nathan

[Daniel Clewley]

Hi Nathan,

If you pass in a list of filenames rather than a single file RIOS will use a list of numpy arrays.

It is documented on the RIOS wiki:

https://bitbucket.org/chchrsc/rios/wiki/applier/listsoffiles

Thanks,

Dan

[Thiago Silva]

Hi guys,

Thank you so much, this is really helpful. Can't believe I missed the stackStats function. And Nathan, I really appreciate you taking the time to write the code on the first reply!

Cheers,

Thiago

To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-support+unsubscribe@googlegroups.com.
To post to this group, send email to rsgislib-support@googlegroups.com.

Visit this group at https://groups.google.com/group/rsgislib-support.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "RSGISLib Support" group.

To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-support+unsubscribe@googlegroups.com.
To post to this group, send email to rsgislib-support@googlegroups.com.