How to calculate multiple equations (a+bX+cY+dZ), where (a, b, and c) are regression coeficients and (X, Y and Z) are rasters, using Pyton?

[LuizVianna]

I need to run 252 equations to calculate freeze probability using lat, lon, and altitude as raster inputs, multiplied by a, b, c and d regression coeficients. The equations are in format a+b*lat+c*lon+d*alt.

The equations are in an Excel spreedsheet (attached) and I need to write them in a python script. The lat, lon and lat raster datasets are in the same folder and the output tif files have a standardized name (also in the Excel file).

Thank you!

[DanPatterson_Retired]

It can

import numpy as np
import arcpy
lat = np.arange(50, 40, -1., 'float')
lon = np.arange(-80, -70, 1., 'float')
lon, lat = np.meshgrid(lon, lat)
tbl = r"C:\GIS\testing_project\testing_project.gdb\equations"
arr = arcpy.da.TableToNumPyArray(tbl, "*")
eqs = arr['EQ']
fns = ["{}{}{}{}".format(*i) for i in arr[['PROB', 'DEC', 'TMIN', 'FORMAT']]]
a = arr['a']
b = arr['b']
c = arr['c']
arrs = [a[i] + b[i]*lat + c[i]*lon for i in range(len(a))]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I have replaced your lon and lat rasters with numpy arrays to demonstrate.  They are formed from a 'meshgrid' using the sample ranges.  They are shown below

The file names are simply put together from your spreadsheet using the ExcelToTable tool in Pro (map has one).

The table was converted to a numpy array (line 7).

You can slice out an equation to check if you want (line 😎

I pieced the bits together to form the file names (line 9)

The constants `a`, `b` `c` are similarly pulled from the array like the other columns

The final results are put together as individual rasters just using a list comprehension (line 13)

Here are some samples

lon  # ---- longitude formed from the meshgrid
 
array([[-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.],
       [-80., -79., -78., -77., -76., -75., -74., -73., -72., -71.]])

lat  # ---- latitude formed from the meshgrid

array([[50., 50., 50., 50., 50., 50., 50., 50., 50., 50.],
       [49., 49., 49., 49., 49., 49., 49., 49., 49., 49.],
       [48., 48., 48., 48., 48., 48., 48., 48., 48., 48.],
       [47., 47., 47., 47., 47., 47., 47., 47., 47., 47.],
       [46., 46., 46., 46., 46., 46., 46., 46., 46., 46.],
       [45., 45., 45., 45., 45., 45., 45., 45., 45., 45.],
       [44., 44., 44., 44., 44., 44., 44., 44., 44., 44.],
       [43., 43., 43., 43., 43., 43., 43., 43., 43., 43.],
       [42., 42., 42., 42., 42., 42., 42., 42., 42., 42.],
       [41., 41., 41., 41., 41., 41., 41., 41., 41., 41.]])


a[20], b[20], c[20]   # ---- 3 cooefficients  for the 21 record
(151.33584980670224, 0.7443302423314689, 2.698509149397158)

arrs[20]   # ---- the 21st array, indexing from 0

array([[-27.3, -24.6, -21.9, -19.2, -16.5, -13.8, -11.1,  -8.4,  -5.7,  -3. ],
       [-28.1, -25.4, -22.7, -20. , -17.3, -14.6, -11.9,  -9.2,  -6.5,  -3.8],
       [-28.8, -26.1, -23.4, -20.7, -18. , -15.3, -12.6,  -9.9,  -7.2,  -4.5],
       [-29.6, -26.9, -24.2, -21.5, -18.8, -16.1, -13.4, -10.7,  -8. ,  -5.3],
       [-30.3, -27.6, -24.9, -22.2, -19.5, -16.8, -14.1, -11.4,  -8.7,  -6. ],
       [-31.1, -28.4, -25.7, -23. , -20.3, -17.6, -14.9, -12.2,  -9.5,  -6.8],
       [-31.8, -29.1, -26.4, -23.7, -21. , -18.3, -15.6, -12.9, -10.2,  -7.5],
       [-32.5, -29.8, -27.1, -24.4, -21.7, -19. , -16.3, -13.6, -11. ,  -8.3],
       [-33.3, -30.6, -27.9, -25.2, -22.5, -19.8, -17.1, -14.4, -11.7,  -9. ],
       [-34. , -31.3, -28.6, -25.9, -23.2, -20.5, -17.8, -15.1, -12.4,  -9.7]])


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

# ---- get the lower left corner from the input lon/lat arrays and make it a point object

LL = arcpy.Point(lon[0,0], lat[-1,0])  # magical slicing

outname = "C:\\Temp\\{}".format(fns[20])

outname
'C:\\Temp\\pgead21n3.tif'

outraster = arcpy.NumPyArrayToRaster(arrs[20], LL, 1., 1.)  # last 2 terms are cell size

outraster.save(outname)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I missed some of the remaining terms, but you get the drift I am sure I have attached the raster as a *.tif

[LuizVianna]

This is what I need:

Follow a data sample.

import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

env.workspace = r"D:\Vianna\Dados\Raster\Clima" #Folder with raster datasets. This line is not working...

#Input raster datasets of the same masked area and same spatial reference
lat = r"D:\Vianna\Dados\Raster\Topografia\latit.tif" #Input lat raster
lon = r"D:\Vianna\Dados\Raster\Topografia\longit.tif" #Input lon raster
alt = r"D:\Vianna\Dados\Raster\Topografia\mde_sc_90.tif" # Input dem raster

tbl = r"D:\Vianna\Dados\Raster\Clima\Raster_Clima.gdb\TB_EQUAC_ETP"

arr = arcpy.da.TableToNumPyArray(tbl, "*")#This line close the ArcGis Pro
eqs = arr['DS_EQUAC']

fns = ["{}{}{}{}".format(*i) for i in arr[['SG_PREFIX', 'SG_SUFIX', 'SG_PERIO', 'SG_FORMA']]]

a = arr['VL_A']
b = arr['VL_B']
c = arr['VL_C']
d = arr['VL_D']

eq = [a[i] + b[i]*lat + c[i]*lon + d[i]*alt for i in range(len(a))] #This is the equation I need

#How to run this equation and write outputs named as "fns"?‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DarrenWiens2]

Try running it in a loop (untested):

arr = arcpy.da.TableToNumPyArray(tbl, "*")

for item in arr:
  fn = "{}{}{}{}".format(
    item['SG_PREFIX'],
    item['SG_SUFIX'],
    item['SG_PERIO'],
    item['SG_FORMA']
  )

  a = item['VL_A']
  b = item['VL_B']
  c = item['VL_C']
  d = item['VL_D']

  eq = a + b*lat + c*lon + d*alt

  ... output here
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[LuizVianna]

The line below closes the  ArcGisPro

arr = arcpy.da.TableToNumPyArray(tbl, "*")

[DanPatterson_Retired]

I made a geodatabase table from the excel file and used it in your line 12

[LuizVianna]

I also did it. And the table is in the file geodatabase. However when I run the line 14, the Python crashes. I´ve tried to run from ArcGis Pro and from prompt.

Did you try to run the script with attached data?

I dont know if the problem is in the script or is the file geodatabase.

[DanPatterson_Retired]

line 6....  tbl = r"C:\GIS\testing_project\testing_project.gdb\equations

I converted your excel to a geodatabase table

I used the code in my example (which you would have to add an extra parameter to I think)

It was for demonstration purposes, but it works

[LuizVianna]

I thank you all, but unfortunatly it still doesn´t working here. I follwed the code step by step, I re-imported the excel table to geodatabase and the line "arr = arcpy.da.TableToNumPyArray(tbl, "*")" still closing Python.

import numpy as np
import arcpy
lat = np.arange(-29.373511, -25.928348, 0.000871310821466946, 'float')
lat
array([-29.373511  , -29.37263969, -29.37176838, ..., -25.93009063,
       -25.92921932, -25.92834801])
lon = np.arange(-53.915588, -48.273851, 0.000871310821466947, 'float')
lon
array([-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
       -48.27559305, -48.27472174])
lon, lat = np.meshgrid(lon, lat)
lat
array([[-29.373511  , -29.373511  , -29.373511  , ..., -29.373511  ,
        -29.373511  , -29.373511  ],
       [-29.37263969, -29.37263969, -29.37263969, ..., -29.37263969,
        -29.37263969, -29.37263969],
       [-29.37176838, -29.37176838, -29.37176838, ..., -29.37176838,
        -29.37176838, -29.37176838],
       ..., 
       [-25.93009063, -25.93009063, -25.93009063, ..., -25.93009063,
        -25.93009063, -25.93009063],
       [-25.92921932, -25.92921932, -25.92921932, ..., -25.92921932,
        -25.92921932, -25.92921932],
       [-25.92834801, -25.92834801, -25.92834801, ..., -25.92834801,
        -25.92834801, -25.92834801]])
lon
array([[-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
        -48.27559305, -48.27472174],
       [-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
        -48.27559305, -48.27472174],
       [-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
        -48.27559305, -48.27472174],
       ..., 
       [-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
        -48.27559305, -48.27472174],
       [-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
        -48.27559305, -48.27472174],
       [-53.915588  , -53.91471669, -53.91384538, ..., -48.27646436,
        -48.27559305, -48.27472174]])
alt = r"D:\Vianna\Dados\Raster\Topografia\mde_sc_90.tif"
tbl = r"D:\Vianna\Dados\Raster\Clima\Raster_Clima.gdb\TB_EQUAC_ETP"
alt
'D:\\Vianna\\Dados\\Raster\\Topografia\\mde_sc_90.tif'
tbl
'D:\\Vianna\\Dados\\Raster\\Clima\\Raster_Clima.gdb\\TB_EQUAC_ETP'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
arr = arcpy.da.TableToNumPyArray(tbl, "*")
eqs = arr['DS_EQUAC']
fns = ["{}{}{}{}".format(*i) for i in arr[['SG_PREFIX', 'SG_SUFIX', 'SG_PERIO', 'SG_FORMA']]]
a = arr['VL_A']
b = arr['VL_B']
c = arr['VL_C']
d = arr['VL_D']
arrs = [a[i] + b[i]*lat + c[i]*lon + d[i]*alt for i in range(len(a))]
LL = arcpy.Point(lon[0,0], lat[-1,0])  # magical slicing
outname = "D:\\Temp\\{}".format(fns[2])
outraster = arcpy.NumPyArrayToRaster(arrs[2], LL, 1., 1.)  # last 2 terms are cell size
outraster.save(outname)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

import arcpy
tbl = r"C:\GIS\testing_project\testing_project.gdb\equations"
arr = arcpy.da.TableToNumPyArray(tbl, "*")
arr.dtype


dtype([('OBJECTID', '<i4'), ('PROB', '<U255'), ('DEC', '<U255'),
 ('TMIN', '<U255'), ('FORMAT', '<U255'), ('a', '<f8'), ('b', '<f8'), ('lat', '<U255'),
 ('c', '<f8'), ('lon', '<U255'), ('d', '<f8'), ('alt', '<U255'), ('EQ', '<U255')])‍‍‍‍‍‍‍‍‍

Just to confirm it works.

Which leaves something wrong with your installation.  You are using ArcGIS Pro's arcpy python (3.6 ish). 

Do you see this as well?

If not, something is wrong with your file and/or the location of the file ie its path... (spaces? on a network, some symlink?)

dir(arcpy.da)
Out[15]: ['Describe', 'Domain', 'Editor', 'ExtendTable', 'FeatureClassToNumPyArray',
 'InsertCursor', 'ListDomains', 'ListFieldConflictFilters', 'ListReplicas',
 'ListSubtypes', 'ListVersions', 'NumPyArrayToFeatureClass', 'NumPyArrayToTable',
 'Replica', 'SearchCursor', 'TableToNumPyArray', 'UpdateCursor', 'Version', 'Walk',
 '__all__', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__',
 '__name__', '__package__', '__spec__', '_internal_eq', '_internal_sd', '_internal_vb']‍‍‍‍‍‍‍

Nothing suggested seems to have helped, so I would contact tech support