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 you are going to have to edit your last code post so it is readable

Code Formatting

plus I don't see where you have converted the input rasters to arrays using RasterToNumPyArray

perhaps show where you have one condition working first before you start with the batch part.  Some may ask for a sample of the three inputs, so I will do it for them.

[XanderBakker]

Below is the code formatted with a bunch of questions. In order to test anything sample data must be provided.

import os
import arcpy
import scipy
import numpy as np

def main():

    arcpy.env.workspace = r"D:\GIS Data\Sixth Creek ET.gdb"
    path1 = r"C:\Workspace\test geowrite\Betasigma"
    path2 = r"C:\Workspace\test geowrite\sigma"
    path3 = r"C:\Workspace\test geowrite\Netrad"

    zGuess= np.array([1,1,1])
    rasters = arcpy.ListRasters("*", "All") # where do you use this?

    lst_files = [f for f in os.listdir(path1) if f.endswith(".tif")]
    lst_fizes = [f for f in os.listdir(path2) if f.endswith(".tif")]
    lst_fixes = [f for f in os.listdir(path2) if f.endswith(".tif")]

    for i in range(0, 46):
        # very risky to do this.
        # does each folder have the same 46 raster?
        # is the order the same so that the correct rasters will be used together?
        path_ras_file = os.path.join(path1, lst_files[i]) 
        path_ras_fize = os.path.join(path1, lst_fizes[i])
        path_ras_fixe = os.path.join(path1, lst_fixes[i])

        # I guess you will have to convert to numpy
        # asuming all have the same dimensiones
        np_file = arcpy.RasterToNumPyArray(arcpy.Raster(path_ras_file))
        np_fize = arcpy.RasterToNumPyArray(arcpy.Raster(path_ras_fize))
        np_fixe = arcpy.RasterToNumPyArray(arcpy.Raster(path_ras_fixe))

        # your function uses a parameter. What is this parameter?
        param ="define parameter!"
        result_function = myFunction(param, np_file, np_fixe, np_fize)

        # are you sure that all the necessary parameters are provided?
        z = scipy.optimize.fsolve(result_function, zGuess)
        
        # do something with the z


def myFunction(z, np_file, np_fixe, np_fize):
    # what is z?
    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 # can this be done with numpy arrays?
    F[1] = np_fixe - E - H - G
    F[2] = np_file * H - E
    return F

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[LarryBiodun]

https://drive.google.com/file/d/1dE32qKVm0V83d4pvmPUD0nWvDmrOaz3w/view?usp=sharing

Hi Xander, thanks for looking at the problem, I have included sample data in my google drive link above (3 each as against 46). A few things are over the top of my head really. I see I don't use the ListRasters command in the script at all, I guess not needed. Yes the rasters are arranged in order i.e the 46 in the same manner in each folder. I don't know about arrays either, just reading about them and the fsolve function. But I am sure the fsolve is able to solve the non linear equations as I tested it earlier using the code below but I replaced the rasters with actual numbers. I am not sure what you mean by parameter in line 32 as I did not add any such in the code below.

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
[  70.96121951  207.89419779   21.14458269]

[RandyBurton]

I'm sure to get in over my head, but it is one way to learn.  So here goes...

I believe  this is a way to pass additional information to fsolve.  It produces the same result as your code example.  I think you are trying to send file, fixe and fize to the function, correct?  These values would come from your raster files.

import numpy as np
from scipy.optimize import fsolve

def myFunction(z, file, fixe, fize):
        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

file = 3.356
fixe = 300
fize = 1.024

zGuess= np.array([1,1,1])
z = fsolve(myFunction, zGuess, args=(file,fixe,fize))
print z

#result: [  70.96121951  207.89419779   21.14458269]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[LarryBiodun]

Hi Randy,

Thanks for taking the shot. I tried this but it returned the same error. Pretty much not wanting to take in arrays or rasters. I have also included some sample raster data above.

[LarryBiodun]

I ran the code you wrote with corrections to the path as shown below and got the error

Runtime error

Traceback (most recent call last):

File "<string>", line 56, in <module>

File "<string>", line 36, in main

File "<string>", line 50, in myFunction

TypeError: only integer arrays with one element can be converted to an index

import os
import arcpy
import scipy
import numpy as np

def main():

    arcpy.env.workspace = r"D:\GIS Data\Sixth Creek ET.gdb"
    path1 = r"C:\Workspace\test geowrite\Betasigma"
    path2 = r"C:\Workspace\test geowrite\sigma"
    path3 = r"C:\Workspace\test geowrite\Netrad"

    zGuess= np.array([1,1,1])
    rasters = arcpy.ListRasters("*", "All") # where do you use this?

    lst_files = [f for f in os.listdir(path1) if f.endswith(".tif")]
    lst_fizes = [f for f in os.listdir(path2) if f.endswith(".tif")]
    lst_fixes = [f for f in os.listdir(path3) if f.endswith(".tif")]

    for i in range(0, 46):
        # very risky to do this.
        # does each folder have the same 46 raster?
        # is the order the same so that the correct rasters will be used together?
        path_ras_file = os.path.join(path1, lst_files[i]) 
        path_ras_fize = os.path.join(path2, lst_fizes[i])
        path_ras_fixe = os.path.join(path3, lst_fixes[i])

        # I guess you will have to convert to numpy
        # asuming all have the same dimensiones
        np_file = arcpy.RasterToNumPyArray(arcpy.Raster(path_ras_file))
        np_fize = arcpy.RasterToNumPyArray(arcpy.Raster(path_ras_fize))
        np_fixe = arcpy.RasterToNumPyArray(arcpy.Raster(path_ras_fixe))

        # your function uses a parameter. What is this parameter?
        param ="define parameter!"
        result_function = myFunction(param, np_file, np_fixe, np_fize)

        # are you sure that all the necessary parameters are provided?
        z = scipy.optimize.fsolve(result_function, zGuess)
        
        # do something with the z


def myFunction(z, np_file, np_fixe, np_fize):
    # what is z?
    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 # can this be done with numpy arrays?
    F[1] = np_fixe - E - H - G
    F[2] = np_file * H - E
    return F

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

it has to be a tuple of 3, notice the comma after the 3

F = np.empty((3,))

[LarryBiodun]

I fixed that and ran the code again and got the same result.

[LarryBiodun]

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)

...

Runtime error

Traceback (most recent call last):

File "<string>", line 28, 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 10, in myFunction

TypeError: unsupported operand type(s) for /: 'str' and 'str'
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[XanderBakker]

Hi Larry Biodun , in order to be more efficient can you work in a IDE and copy the code without the indent points and error messages to have only the code preserving the indentation? It would require quite some editing to be able to run the code and the idea is when you insert code that it can be copied and executed, which is not the case right now. Just have a look at the code I inserted earlier.