topic Re: Logarithmic Averaging / Geometric Mean in ArcGIS Spatial Analyst Questions

Logarithmic Averaging / Geometric Mean

KoreenMillard — Mon, 04 Oct 2010 17:33:54 GMT

Hello, I want to calculate the logarithmic average of the values of points as a raster.

Has anyone ever tried to do this?

-k

FYI:

Geometric mean explained by wikipedia: "The geometric mean, in mathematics, is a type of mean or average, which indicates the central tendency or typical value of a set of numbers. It is similar to the arithmetic mean, which is what most people think of with the word "average", except that the numbers are multiplied and then the nth root (where n is the count of numbers in the set) of the resulting product is taken."

http://en.wikipedia.org/wiki/Geometric_mean

Re: Logarithmic Averaging / Geometric Mean

KoreenMillard — Fri, 08 Oct 2010 14:28:22 GMT

Thanks very much!

Re: Logarithmic Averaging / Geometric Mean

MichaelAugust — Fri, 31 Oct 2014 18:16:10 GMT

Hi Koreen - did you figure out how to calculate the geometric mean with the raster calculator? Thanks!

Re: Logarithmic Averaging / Geometric Mean

curtvprice — Thu, 06 Nov 2014 07:26:49 GMT

The way you'd do this using the 10x Raster Calculator tool would be:

("grid1" * "grid2" * "grid3" ) ** (1 / 3.)

Don't forget to include the "." after the number of rasters (3), as Python division is integer if both inputs are integer.

>>> 1 / 3
0
>>> 1 / 3.
0.3333333333333333

Re: Logarithmic Averaging / Geometric Mean

MichaelAugust — Thu, 06 Nov 2014 22:13:29 GMT

Hi Curtis - I meant calculate the geometric mean of all the individual cells contained in one raster, like this:

with each of those being an individual cell of my raster. I think we realized that doing this to 1.2M cells in a raster simply implodes in Python, we couldn't get the script to do more than 50 or so values at one time. Can someone out there with more programming experience explain why, or if there's a better way to do geometric mean in Python or ArcGIS?
thanks!

Re: Logarithmic Averaging / Geometric Mean

ShaunWalbridge — Thu, 06 Nov 2014 22:59:39 GMT

Not a spatial problem per-se, and the SciPy library does a pretty good job for problems like this. It'll be included in the 10.3 release, but can be downloaded and installed separately in 10.x. With that in place, convert your raster to a NumPy array, then use the statistics built in to SciPy to do the calculation:

import scipy.stats
rast_path = 'C:/my_input_raster.tif')
raster_as_numpy_array = arcpy.RasterToNumPyArray(rast_path)
raster_geometric_mean = scipy.stats.stats.gmean(
    raster_as_numpy_array, axis=None)

Which for a raster of 10M values takes a few seconds for me to run.

Hope that helps,

Shaun

Edit: Michael August I've updated my answer to include 'axis=None' which will compute the geometric mean for the whole matrix (instead of along one axis). With this change, it should work without any further steps.

Re: Logarithmic Averaging / Geometric Mean

DanPatterson_Retired — Sat, 11 Dec 2021 20:57:13 GMT

This is an alternate example of Curtis's example, not specifically what the OP wants. ...In the example below I created 3 numpy arrays representing your example, assuming everything under the root sign is a grid.


>>> import numpy as np


>>> base_form = np.arange(0,9,dtype='float64')


>>> base_form.reshape((3,3))


array([[ 0.,  1.,  2.],


      [ 3.,  4.,  5.],


      [ 6.,  7.,  8.]])


>>> a = ((np.ones_like(base_form))*2).reshape((3,3))


>>> b = ((np.ones_like(base_form))*4).reshape((3,3))


>>> c = ((np.ones_like(base_form))*8).reshape((3,3))


>>>


>>> d = a*b*c


>>> d


array([[ 64.,  64.,  64.],


      [ 64.,  64.,  64.],


      [ 64.,  64.,  64.]])


>>> root = ((np.ones_like(base_form))*(1.0/3.0)).reshape((3,3))


>>> root


array([[ 0.33333333,  0.33333333,  0.33333333],


      [ 0.33333333,  0.33333333,  0.33333333],


      [ 0.33333333,  0.33333333,  0.33333333]])


>>> gm = d**root


>>> gm


array([[ 4.,  4.,  4.],


      [ 4.,  4.,  4.],


      [ 4.,  4.,  4.]])


>>>

it might be worth considering since you can use Raster to Numpy Array and the reverse to go between numpy and Arc

Now on to the simpler question for the gm of a single raster. You need the cumulative product of each cell, getting the last value in the row column which you do with a np.flatten function if it is a n*m array. For a simple array. the syntax would be


>>> # different example
>>> a = (np.arange(1,10)).reshape((3,3))
>>> a
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> cp = np.cumproduct(a)
>>> cp
array([     1,      2,      6,     24,    120,    720,   5040,  40320,
       362880])
>>> N = np.size(cp)
>>> N
9
>>> gm = (cp[-1])**(1.0/N)
>>> gm
4.1471662743969127
>>>


Which is the same as excel...for your example,


on just those numbers....


>>> a = np.array([2.,4.,8.])


>>> cp = np.cumprod(a)


>>> cp


array([  2.,   8.,  64.])


>>> gm = cp[-1]**(1.0/np.size(cp))


>>> gm


3.9999999999999996


>>>

You can also use the alternate formulation to calculate the geometric mean for large numbers and/or arrays to prevent numeric over/under flows. For example, a case of the 3x3 array as before and 1,000,000 int32 integers between 1 and 255


>>> import numpy as np


>>> import random


>>> #************************ Simple example using log10  **************************


>>> a = (np.arange(1,10)).reshape((3,3))


>>> loga = np.log10(a)


>>> loga


array([[ 0.        ,  0.30103   ,  0.47712125],


       [ 0.60205999,  0.69897   ,  0.77815125],


       [ 0.84509804,  0.90308999,  0.95424251]])


>>> cum_sum = np.cumsum(loga)


>>> N = np.size(a)


>>> N


9


>>> gm = 10.0**((cum_sum[-1])/N)


>>> gm


4.1471662743969135


>>> import numpy as np


>>> import random


>>> #***************************  large array * **************


>>> big = [ random.randint(1,255) for i in xrange(1000000)]


>>> a = np.array(big,dtype='int64')


>>> loga = np.log10(a)


>>> loga


array([ 2.32221929,  1.95904139,  2.38201704, ...,  2.19865709,


        1.68124124,  2.25285303])


>>> N = len(big)


>>> N


1000000


>>> cum_sum = np.cumsum(loga)


>>> cum_sum


array([  2.32221929e+00,   4.28126069e+00,   6.66327773e+00, ...,


         1.97970102e+06,   1.97970270e+06,   1.97970495e+06])


>>> gm = 10.0**((cum_sum[-1])/N)


>>> gm


95.434401572885349


>>>

An no other packages need to be installed

Re: Logarithmic Averaging / Geometric Mean

curtvprice — Fri, 07 Nov 2014 05:36:48 GMT

Once upon a time there was a tool in ArcInfo Workstation GRID called DOCELL that would operate through all the cells of a raster. It was cool:

Grid: g = scalar(1)
Grid: n = scalar(0)
Grid: DOCELL
:: g *= ingrid
:: n += 1
::  END
Grid: g ^ (1 / float(n))

Talk about a blast from the past.

Shaun Walbridge‌ 's suggestion like a good way to do it nowadays. But if you can't get scipy installed (a bit of a lift actually) you could do something similar by converting the raster to a numpy array with RasterToNumpyArray() and doing the calculation there. That's what Dan did.

I just had a another thought - if the input is an integer grid with a VAT, something like this would work:

sum = 0
ncells = 0
with arcpy.da.SearchCursor("ingrid", ["VALUE", "COUNT"]) as rows:
    for row in rows:
        for cc in range(row[1]): 
            sum += sum * row[0] # running geometric sum of all cells in zone
        ncells += row[1]        # running cell count
gmean = sum ** (1 / float(ncells))

It's possible you may run into nasty numerical problems if the grid is very big as the sum gets huge, but this avoids converting the raster to numpy.

Re: Logarithmic Averaging / Geometric Mean

MichaelAugust — Fri, 07 Nov 2014 18:47:05 GMT

>>> raster_as_numpy_array = arcpy.RasterToNumPyArray('C:\Users\xxxxx\Downloads\q47121g52be.tif\q47121g52be.tif')
>>> raster_geometric_mean = scipy.stats.mstats.gmean(raster_as_numpy_array)
>>> raster_geometric_mean
array([ 0.,  0.,  0., ...,  0.,  0.,  0.], dtype=float16)

So I tried the above on integer and float rasters and am not sure what I did wrong there...how do you get to the gm as a single reported value? If you can't tell I am fairly new at this point to Python -thanks you guys for all your help!

Re: Logarithmic Averaging / Geometric Mean

DanPatterson_Retired — Fri, 07 Nov 2014 18:55:26 GMT

once you have the array...check my examples in this thread

Re: Logarithmic Averaging / Geometric Mean

MichaelAugust — Fri, 07 Nov 2014 19:01:17 GMT

So finish it with this part:

>>> cp = np.cumproduct(a) 

>>> cp 
array([    1,      2,      6,    24,    120,    720,  5040,  40320,  
362880])  
>>> N = np.size(cp) 
>>> N 
9  
>>> gm = (cp[-1])**(1.0/N)  
>>> gm 
4.1471662743969127

Re: Logarithmic Averaging / Geometric Mean

DanPatterson_Retired — Fri, 07 Nov 2014 19:07:15 GMT

In my thread, I outlined two options...one using the cumulative product approach...and the other..the cumulative sum of the logs... the cumulative product approach will fail if you have a large array and/or the numbers in the array contain large numbers due to numeric overflow. The log approach suffers less from this. The only difference is the way in which the final GM is determined. I documented this for teaching purposes and to provide some context on GM use and the calculation pitfalls. I am pretty sure that spreadsheet geometric mean functions don't use the cumulative product approach but rollout the log summation method behind the scenes..... So read my thread carefully...

Re: Logarithmic Averaging / Geometric Mean

MichaelAugust — Tue, 18 Nov 2014 20:36:34 GMT

Great - thanks Shaun!