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?

5687
35
Jump to solution
09-03-2018 01:23 PM
LuizVianna
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
LuizVianna
Occasional Contributor

Thank you all!

The final code is:

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

wksp = "D:\Vianna\Dados\Raster\Clima\ETP"
arcpy.env.workspace = wksp

lat = arcpy.Raster("D:\Vianna\Dados\Raster\Topografia\latit.tif")
lon = arcpy.Raster("D:\Vianna\Dados\Raster\Topografia\longit.tif")
alt = arcpy.Raster("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']
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])

View solution in original post

35 Replies
DanPatterson_Retired
MVP Emeritus

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

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"?‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
DarrenWiens2
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
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
LuizVianna
Occasional Contributor

The line below closes the  ArcGisPro

arr = arcpy.da.TableToNumPyArray(tbl, "*")
0 Kudos
DanPatterson_Retired
MVP Emeritus

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

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

0 Kudos
DanPatterson_Retired
MVP Emeritus

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

0 Kudos
LuizVianna
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 = 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)
0 Kudos
DanPatterson_Retired
MVP Emeritus
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