topic Re: Solving 3 non linear algebraic equations in python in Python Questions

Solving 3 non linear algebraic equations in python

LarryBiodun — Wed, 29 Nov 2017 04:44:55 GMT

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.

Re: Solving 3 non linear algebraic equations in python

DanPatterson_Retired — Wed, 29 Nov 2017 11:12:20 GMT

If you can solve it with SciPy and the input takes 2D arrays then you can convert the raster to an array using

RasterToNumPyArray

the only concern is to ensure that the extent and cell sizes are the same. If the data types are the same that is good, but the array can be promoted or demoted as necessary.

If the inputs require a 1D array, then you can 'ravel()' the 2D array.

Perhaps some more information would facilitate devising a workflow to permit the integration

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Wed, 29 Nov 2017 12:57:30 GMT

Hi Dan,

Thanks for your response.

While I thought it should be fairly easy with scipy...I was wrong. I wrote this basic part of the code however it didn't quite work with all the errors below and I am struggling to figure it out. As u can guess I am a newbie.

My struggles are now that I cant even get the fsolve to work. I found out scipy is included in ArcGIS 10.4.1 (which is my current version), so I imported it and tried my hands on it. As a sample instead of trying to work with rasters, I used single values for the raster component...with no success. The values 1.024 in F(0), 300 in F(1) and 3.356 in F(0) and F(2) are daily rasters which I want the code to iterate through and pick the set for each day. There are 46 each in different folders. I hope the information below can help explain the problems a bit more clearly. Thanks.

import scipy

import numpy as np

def myFunction(z):
...     E = z[0]
...     G = z[1]
...     H = z[2]
...     F = 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)
Traceback (most recent call last):

File "<ipython-input-29-eb6f86348f14>", line 1, in <module>
z = fsolve(myFunction,zGuess)

File "C:\Program Files\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 146, in fsolve
res = _root_hybr(func, x0, args, jac=fprime, **options)

File "C:\Program Files\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 212, in _root_hybr
shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))

File "C:\Program Files\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 26, in _check_func
res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))

File "<ipython-input-27-7e512d399a9a>", line 5, in myFunction
F = empty((3))

NameError: name 'empty' is not defined

Re: Solving 3 non linear algebraic equations in python

DanPatterson_Retired — Wed, 29 Nov 2017 14:23:13 GMT

immediate issues... I will deal with the rest after classes

empty... can't just doe that, it has to be np.empty

when you import numpy as np... as other functions should be preceeded by np... don't do an import all either to save two letterd

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Wed, 29 Nov 2017 14:43:22 GMT

Right! Makes sense. Noted. Thanks I'll wait for your response after your class.

Cheers

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Wed, 29 Nov 2017 23:17:25 GMT

Hi Dan, I fixed the numpy thing as you said. I fiddled around with it a little and it appears to be working now, however, I still need help with coding it to take those rasters in sequentially and also to write each of the E, G and H out as individual raster for each time step it iterates through. See below the code that worked. Look forward to your response.

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]

Re: Solving 3 non linear algebraic equations in python

DanPatterson_Retired — Fri, 10 Dec 2021 20:23:28 GMT

just got back fsolve is in the scipy optimize minpack

https://github.com/scipy/scipy/blob/05028ff66eadeee32b33ac2f994c009093355534/scipy/optimize/minpack.py

notice how I imported fsolve on line 1.

I cut out some of the 'dir' information on line 4

The full help is on line 9 onwards.

So on first glance,it you have the 3 parameters then the tuple of 3 seems reasonable. Interpreting it I suspect would be in the order that they were given.

If will look further soon, but this should keep you up to date

from scipy.optimize import fsolve

dir(fsolve)
Out[6]: ['__annotations__', '__call__', '__class__...._name__', '__ne__', '__new__', '__qualname__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__']


help(fsolve)

Help on function fsolve in module scipy.optimize.minpack:

fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)
    Find the roots of a function.
    
    Return the roots of the (non-linear) equations defined by
    ``func(x) = 0`` given a starting estimate.
    
    Parameters
    ----------
    func : callable ``f(x, *args)``
        A function that takes at least one (possibly vector) argument.
    x0 : ndarray
        The starting estimate for the roots of ``func(x) = 0``.
    args : tuple, optional
        Any extra arguments to `func`.
    fprime : callable(x), optional
        A function to compute the Jacobian of `func` with derivatives
        across the rows. By default, the Jacobian will be estimated.
    full_output : bool, optional
        If True, return optional outputs.
    col_deriv : bool, optional
        Specify whether the Jacobian function computes derivatives down
        the columns (faster, because there is no transpose operation).
    xtol : float, optional
        The calculation will terminate if the relative error between two
        consecutive iterates is at most `xtol`.
    maxfev : int, optional
        The maximum number of calls to the function. If zero, then
        ``100*(N+1)`` is the maximum where N is the number of elements
        in `x0`.
    band : tuple, optional
        If set to a two-sequence containing the number of sub- and
        super-diagonals within the band of the Jacobi matrix, the
        Jacobi matrix is considered banded (only for ``fprime=None``).
    epsfcn : float, optional
        A suitable step length for the forward-difference
        approximation of the Jacobian (for ``fprime=None``). If
        `epsfcn` is less than the machine precision, it is assumed
        that the relative errors in the functions are of the order of
        the machine precision.
    factor : float, optional
        A parameter determining the initial step bound
        (``factor * || diag * x||``).  Should be in the interval
        ``(0.1, 100)``.
    diag : sequence, optional
        N positive entries that serve as a scale factors for the
        variables.
    
    Returns
    -------
    x : ndarray
        The solution (or the result of the last iteration for
        an unsuccessful call).
    infodict : dict
        A dictionary of optional outputs with the keys:
    
        ``nfev``
            number of function calls
        ``njev``
            number of Jacobian calls
        ``fvec``
            function evaluated at the output
        ``fjac``
            the orthogonal matrix, q, produced by the QR
            factorization of the final approximate Jacobian
            matrix, stored column wise
        ``r``
            upper triangular matrix produced by QR factorization
            of the same matrix
        ``qtf``
            the vector ``(transpose(q) * fvec)``
    
    ier : int
        An integer flag.  Set to 1 if a solution was found, otherwise refer
        to `mesg` for more information.
    mesg : str
        If no solution is found, `mesg` details the cause of failure.
    
    See Also
    --------
    root : Interface to root finding algorithms for multivariate
    functions. See the 'hybr' `method` in particular.
    
    Notes
    -----
    ``fsolve`` is a wrapper around MINPACK's hybrd and hybrj algorithms.

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Wed, 29 Nov 2017 23:49:45 GMT

Yes Dan thanks. I forgot to include the line that I used to import the fsolve function in scipy. I also agree that the results are written in the order of input. From here I guess it's how to iterate through through the folders to pick the inputs and also write out the outputs as rasters that will seal the work!

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Thu, 30 Nov 2017 06:14:18 GMT

After getting the code to run using just numbers to replace the rasters, I tried this code but it didn't quite work. Please help. This is doing my head in. I tried this code which I basically tried to cut out a chunk of a code that once worked in iteration. Please help. I am running out of time for a dealine

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

... fize = path2+fizes

... fixe = path3+fixes

... 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'

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Thu, 30 Nov 2017 06:23:56 GMT

Xander Bakker,@Dan Patterson, Randy burton, @kevin Dunlop, @Joshua Bixby After getting the code to run using just numbers to replace the rasters, I tried this code but it didn't quite work. Please help. This is doing my head in. I tried this code which I basically tried to cut out a chunk of a code that once worked in iteration. Please help. I am running out of time for a dealine

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

... fize = path2+fizes

... fixe = path3+fixes

... 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'

Re: Solving 3 non linear algebraic equations in python

DanPatterson_Retired — Thu, 30 Nov 2017 07:04:11 GMT

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.

Re: Solving 3 non linear algebraic equations in python

XanderBakker — Fri, 10 Dec 2021 20:23:32 GMT

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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Thu, 30 Nov 2017 15:24:38 GMT

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.

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

print z
[ 70.96121951 207.89419779 21.14458269]

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Fri, 10 Dec 2021 20:23:34 GMT

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'
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Solving 3 non linear algebraic equations in python

XanderBakker — Thu, 30 Nov 2017 16:16:20 GMT

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.

Re: Solving 3 non linear algebraic equations in python

DanPatterson_Retired — Thu, 30 Nov 2017 16:21:56 GMT

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

Re: Solving 3 non linear algebraic equations in python

XanderBakker — Thu, 30 Nov 2017 16:49:50 GMT

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.

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Fri, 10 Dec 2021 20:23:37 GMT

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‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Fri, 10 Dec 2021 20:23:40 GMT

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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Solving 3 non linear algebraic equations in python

LarryBiodun — Fri, 10 Dec 2021 20:23:43 GMT

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()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍