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
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 <SPAN class="o" style="font-size: 11.1999998092651px;">=</SPAN> <SPAN class="n" style="font-size: 11.1999998092651px;">arcpy</SPAN><SPAN class="o" style="font-size: 11.1999998092651px;">.</SPAN><SPAN class="n" style="font-size: 11.1999998092651px;">RasterToNumPyArray(rast_path</SPAN><SPAN class="p" style="font-size: 11.1999998092651px;">)</SPAN> raster_geometric_mean = scipy.stats.stats.gmean( raster_as_numpy_array, axis=None)
<SPAN class="o" style="font-size: 11.1999998092651px;">=</SPAN> <SPAN class="n" style="font-size: 11.1999998092651px;">arcpy</SPAN><SPAN class="o" style="font-size: 11.1999998092651px;">.</SPAN><SPAN class="n" style="font-size: 11.1999998092651px;">RasterToNumPyArray(rast_path</SPAN><SPAN class="p" style="font-size: 11.1999998092651px;">)</SPAN>
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.
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!
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
Hi Koreen - did you figure out how to calculate the geometric mean with the raster calculator? Thanks!
Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.
