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?

3169
35
09-03-2018 01:23 PM
Occasional Contributor

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!

1 Solution

Accepted Solutions
Occasional Contributor

Thank you all!

The final code is:

``````import numpy as np
import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True

arcpy.env.workspace = wksp

arr = arcpy.da.TableToNumPyArray(tbl, "*")
eqs = arr['DS_EQUAC']
filenames = ["{}{}{}{}".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']

for i in range(len(filenames)):
print(filenames[i])
eq = a[i] + b[i] * lat + c[i] * lon + d[i] * alt
eq.save(filenames[i])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````
35 Replies
MVP Esteemed Contributor

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

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

Occasional Contributor

This is what I need:

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

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"?‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````
MVP Honored Contributor

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
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````
Occasional Contributor

The line below closes the  ArcGisPro

``arr = arcpy.da.TableToNumPyArray(tbl, "*")``
MVP Esteemed Contributor

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

Occasional Contributor

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.

MVP Esteemed Contributor

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

Occasional Contributor

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
tbl
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)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````
MVP Esteemed Contributor
``````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