Solving 3 non linear algebraic equations in python

[LarryBiodun]

How can I solve a non-linear algebraic equation in ArcGIS python over multiple rasters. Let me Rephrase. I want to solve the following 3 non linear equations ,    and for 46 8 day time steps.  I have 46 rasters each for an 8 day period for Β(σ) , and σ, where I need to take input values from per time step. Is/Io is a constant. The three unknowns I have are E, G and H. The solution should ideally produce 46 rasters each for the 3 unknowns. I don't know if this can be done in ArcGIS though. I know I can solve this equations easily in SciPy using the Fsolve function but I don't know how to get this working with rasters. Any code that could help will be greatly appreciated.

[DanPatterson_Retired]

Larry appears to be running it interactively rather than running a script from an IDE

[XanderBakker]

I noticed that but if you simply copy and paste the code, it will not be directly usable in an IDE (or anywhere else), hence the suggestion to use an IDE. He could also right click and save the code and share the relevant part to avoid all the additional points to be included and to loose the indentation. 

[LarryBiodun]

I think I understand what you mean. I'll post both codes I've worked with in the format you requested. First the little code that ran the equation and printed values for E,H and G.

import arcpy
import scipy
import numpy as np
def myFunction(z):
        E = z[0]
        G = z[1]
        H = z[2]
        F = np.empty((3))
        F[0] = (3.356/1.024)*3*H*pow(abs(H),-(1/6))-G
        F[1] = 300-E-H-G
        F[2] = 3.356*H-E
        return F

zGuess= np.array([1,1,1])
z = fsolve(myFunction,zGuess)

print z‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[LarryBiodun]

Sorry for the late response, it's still in the very early hours of the morning in my part of the world. I guess this is what you requested. The code below was giving me error when I ran it.

import os

import arcpy

import scipy

import numpy as np

def myFunction(z):

    E = z[0]

    G = z[1]

    H = z[2]

    F = np.empty((3))

    F[0] = (file/fize)*3*H*pow(abs(H),-(1/6))-G

    F[1] = fixe-E-H-G

    F[2] = file*H-E

    return F



zGuess= np.array([1,1,1])

arcpy.env.workspace ="D:/GIS Data/Sixth Creek ET.gdb"

rasters = arcpy.ListRasters("*", "All")

path = "C:\\Workspace\\test geowrite\\Betasigma\\"

path2 = "C:\\Workspace\\test geowrite\\sigma\\"

path3 = "C:\\Workspace\\test geowrite\\Netrad\\"

files= [f for f in os.listdir(path) if f.endswith(".tif")]

fizes= [f for f in os.listdir(path2) if f.endswith(".tif")]

fixes= [f for f in os.listdir(path2) if f.endswith(".tif")]

for i in range(0,46):

    file = path+files[i]

    fize = path2+fizes[i]

    fixe = path3+fixes[i]

    z = fsolve(myFunction,zGuess)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[LarryBiodun]

Xander BakkerDan Patterson think I am getting closer to concluding that the Fsolve function may not be capable of doing what I need it to do.

I ran my code for a single set of rasters below and got the error below. It appears that just as Xander Bakker suggested in line 50 of his code, maybe this cant be done with numpy arrays.  Any other ideas of how to get this solved?

Runtime error 
Traceback (most recent call last):
File "<string>", line 29, in <module>
File "C:\Program Files (x86)\python27\ArcGIS10.4\lib\site-packages\scipy\optimize\minpack.py", line 140, in fsolve
res = _root_hybr(func, x0, args, jac=fprime, **options)
File "C:\Program Files (x86)\python27\ArcGIS10.4\lib\site-packages\scipy\optimize\minpack.py", line 197, in _root_hybr
shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
File "C:\Program Files (x86)\python27\ArcGIS10.4\lib\site-packages\scipy\optimize\minpack.py", line 20, in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
File "<string>", line 23, in myFunction
ValueError: setting an array element with a sequence.

import os
import arcpy
import scipy
import numpy as np
 
def myFunction(z):
     inRas1 = arcpy.Raster("C:\\Workspace\\test geowrite\\Betasigma\\Rad_meantiff_001.tif")
     lowerleft = arcpy.Point(inRas1.extent.XMin,inRas1.extent.YMin)
     cellSize = inRas1.meanCellWidth
     np_file = arcpy.RasterToNumPyArray(inRas1,nodata_to_value=0)
     inRas2 = arcpy.Raster("C:\\Workspace\\test geowrite\\sigma\\Rad_meantiff_001.tif")
     lowerleft = arcpy.Point(inRas2.extent.XMin,inRas2.extent.YMin)
     cellSize = inRas2.meanCellWidth
     np_fize = arcpy.RasterToNumPyArray(inRas2,nodata_to_value=0)
     inRas3 = arcpy.Raster("C:\\Workspace\\test geowrite\\Netrad\\Rad_meantiff_001.tif")
     lowerleft = arcpy.Point(inRas3.extent.XMin,inRas3.extent.YMin)
     cellSize = inRas3.meanCellWidth
     np_fixe = arcpy.RasterToNumPyArray(inRas3,nodata_to_value=0)     
     E = z[0]
     G = z[1]
     H = z[2]
     F=np.empty((3,))
     F[0] = (np_file/np_fize)*3*H*pow(abs(H),-(1/6))-G
     F[1] = np_fixe-E-H-G
     F[2] = np_file*H-E
     return F
      
zGuess= np.array([1,1,1])
z = fsolve(myFunction,zGuess)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

my only thought at the moment is that you are using 2D arrays as input.  I suspect that fsolve is looking for a single value for each of the inputs... effectively burrowing downwards through the stack of inputs and it only solves for each combination of the inputs. 

What you are trying to do is get a value for each of the inputs at each location.  This would be analogous to doing the root solve for each location based on the values at each location.  It doesn't seem to be able to handle the 2D inputs but requires a list/tuple of the values for zguess at each location

[LarryBiodun]

I think I get your point. However I have no idea on how to possibly attempt to extract one variable at a time. Could you possibly write something that I can insert into the code to do the suggestion you made?

[DanPatterson_Retired]

visually what I am talking about and how the values need to be generated... they do not represent the actual data you are using but represent 3 arrays with values and how do you generate the inputs from them.

# ---- 3 arrays, a, b, c ----
a 
array([[1, 4, 1],
       [1, 3, 1],
       [2, 4, 2]])

b
array([[1, 3, 1],
       [2, 4, 2],
       [1, 4, 1]])

c
array([[8, 7, 6],
       [6, 8, 8],
       [6, 8, 7]])

# ---- flatten them out so you can have a look ------
a0 = a.ravel(); b0 = b.ravel(); c0 = c.ravel()

a0, b0, c0
(array([1, 4, 1, 1, 3, 1, 2, 4, 2]),
 array([1, 3, 1, 2, 4, 2, 1, 4, 1]),
 array([8, 7, 6, 6, 8, 8, 6, 8, 7]))

# ---- now look down each 'column' aka location to get the values----
d = [[a0[i], b0[i], c0[i]] for i in range(0, 9)]

d 
[[1, 1, 8],
 [4, 3, 7],
 [1, 1, 6],
 [1, 2, 6],
 [3, 4, 8],
 [1, 2, 8],
 [2, 1, 6],
 [4, 4, 8],
 [2, 1, 7]]
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Which basically means the inputs can be reshaped

d0 = np.asarray(d)

d0.reshape(3,3,3)
 
array([[[1, 1, 8],
        [4, 3, 7],
        [1, 1, 6]],

       [[1, 2, 6],
        [3, 4, 8],
        [1, 2, 8]],

       [[2, 1, 6],
        [4, 4, 8],
        [2, 1, 7]]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

So the original inputs from array a, b, c can be reshuffled in a variety of ways, but the formulation of 'd' is probably easier to visualize how to get 3 values for each location ... aka combination, of a, b, and c.

In a practical sense, you would need to find out the plausible combinations of your inputs to produce a 'lookup' table

[LarryBiodun]

Please this is still an on going challenge. If you have any other suggestions. It will be very much appreciated. The flattening of the arrays did not work. 

[DanPatterson_Retired]

Still think the lookup table with the 3 values as input and one value as output.

You would have to produce that table using 'Combine'

You need to indicate how many different classes do you have for each variable?  That will tell you the combinations possible.

By all appearance there is no way to solve the whole array at once... an array of combinations will give you the output value for each triplet.  The table of triplets and result gives you the lookup table which can be used to reclass the combined array (using Combine) which is the combinations of the 3 variables at each location