calculations on rasters stored in memory

[MattiasVan_Opstal]

Hi everyone,

I want to write a program, in which I do calculations on several rasters.  In the program, I use some rasters that I calculated earlier in (the same) program. The problem is now that this doesn't seem to work.
The program crashes and in my geodatabase a 'timesras, ifthe_ras,... are found after the crash.
Here follows a small piece of code with the same logic as what I do in my program, what do I do wrong here? I hope I don't have to save each raster (here: a and b) before I can use them in a new calculation? Are there other sollutions?



Kind Regards,
Mattias

import arcpy
from arcpy import *
from arcpy.sa import *

env.workspace = "C:/..."
arcpy.CheckOutExtension("Spatial")

InRas1 = Raster('C:/...')
InRas2 = Raster('C:/...')

a = InRas1*InRas2
b = InRas1/InRas2

OutRas = a*b
OutRas.save('C:/...')

[Luke_Pinner]

That piece of code ran fine for me when when I changed the file paths to actual rasters.  Make sure your paths don't have any spaces in them.

[MattiasVan_Opstal]

Thank you for the reply, but I still have the problem, is it possible that the problem only occurs on larger rasters/with more calculations?
What I also have (for this example) is that the intermediate rasters are found in my geodatabase after I ran the code (times_ras and divid_ras), any way to avoid this?

[JeffreySwain]

Also to echo what Luke mentioned, be sure you are using the proper pathing.  In your example, you are not using the r".." nor the \\ necessary for python paths.  I am not sure if that is causing the issue or not.  Here is the help documentation that describes this.  But if the paths are correct, then the only other issue maybe how much overlap and the coordinate system of the input rasters.  I recommend making them the same coordinate system and like all other raster processing, having a common extent and not trying to multiply/divide NULL values is important.  I would not think that the size of the rasters would cause an issue, unless they are so huge that you could not process them normally if you broke down the process step by step.  Can you indicate what error message you are receiving? Or is it a direct crash?

I would recommend changing the working environment if you are getting intermediate results back in the geodatabase. There is information on setting the environment in python as well.

[DanPatterson_Retired]

Jeff a note of correction
Python uses \\ in paths or a single / or the raw notation if a single backslash is being used.

[MattiasVan_Opstal]

I am using the correct pathing. I also don't think it has something to do with the size of the rasters, the problem stays the same when I work with small rasters.
What I found out now, is that when I step trough with the debugger in PythonWin, I don't get an errormessage (I can check the locations of the intermediate rasters in the command window while stepping trough, all goes fine). If I do it like this, the output raster will be saved at the end of the program without a problem. When
I run the program without debugging, it crashes (virtual function call) at the moment I try to save the raster (before all goes fine)

In the program, I do around 50 calculations (most of them containing complex Con() statements) on different rasters that are placed in 2 different classes

[MattiasVan_Opstal]

Maybe it is easier if I give a piece of the code, here I added the smallest piece of code that gives a problem.
If I run this, the program will crash, giving a runtime error, Pure virtual function call. However, if I uncomment the intermediate saves,  there is no problem anymore and the progam will work fine. I get the feeling it has something to do with the memory use of the Con statements... Sometimes I nest 3 Con-statements into each other, within the Con statements there are functions containing again 3 Con-statements themself...

If someone is willing to test this piece of code, for the first 6 rasters I use random rasters (values between 0 and 1) of 10 by 10 cells. The ZeroRaster is a 10 by 10 raster containing only zero's.

Thank you!
Mattias
[PHP]
# Import moduls

import arcpy
from arcpy import *
from arcpy.sa import *

# Set environment settings

env.workspace = "C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb"

# Check out the ArcGIS Spatial Analyst extension license

arcpy.CheckOutExtension("Spatial")

# Set variables

K_sat_L1 = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/K_sat_L1xx')
theta_C1 = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/theta_C1xx')
theta_fc_L1 = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/theta_fc_L1xx')
theta_sat_L1 = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/theta_sat_L1xx')
thickness_C1 = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/thickness_C1xx')
tau_L1 = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/tau_L1xx')
ZeroRaster = Raster('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/ZeroRasterxx')


                    #----soil compartments----#

# class for soil compartments

class SoilCompartment:
    def __init__(self, theta, thickness,theta_sat,theta_fc,K_sat,tau):
       
        self.theta = theta              
        self.thickness = thickness
        self.theta_sat = theta_sat
        self.theta_fc = theta_fc
        self.K_sat = K_sat
        self.tau = tau

# create array of SoilCompartment objects
       
Compartment = []
Compartment.append(SoilCompartment(theta_C1,thickness_C1,theta_sat_L1,theta_fc_L1,K_sat_L1,tau_L1)) # first compartment


#given theta -> calculate the corresponding delta theta/ delta t


def calculate_delta_theta(theta, theta_sat, theta_fc, tau):
# theta > theta_sat:
    delta_theta = Con(theta > theta_sat, theta_sat,
# theta < theta_fc:
    Con(theta <= theta_fc, 0,
# theta_fc < theta < theta_sat:
    tau*(theta_sat - theta_fc)*(Exp(theta-theta_fc)-1)/(Exp(theta_sat - theta_fc)-1)))
    return delta_theta

#given delta theta/ delta t -> calculate the corresponding theta

def calculate_theta(delta_theta, theta_sat, theta_fc, tau):
#delta_theta <= 0
    theta = Con(delta_theta <= 0,theta_fc,
#tau > 0
    Con(tau > 0, theta_fc + Log2(1 + delta_theta*(Exp(theta_sat - theta_fc) - 1))/(tau*(theta_sat - theta_fc)),
#else, to stop draining:
    theta_sat + 0.1))
    return theta


                       #----calculation Procedure for drainage----#

drainsum = ZeroRaster                  
delta_theta = ZeroRaster
drain_comp = ZeroRaster
drainmax = ZeroRaster
theta_x = ZeroRaster
excess = ZeroRaster
pre_thick = ZeroRaster

#calculate drainage

#Drainsum cannot exceed the K_sat of the soil layer

def CheckExcess(drainsum, excess, K_sat):
    excess = Con(drainsum > K_sat, excess + drainsum - K_sat,
    excess)
    return excess
   
def CheckDrainsum(drainsum, K_sat):
    drainsum = Con(drainsum > K_sat, K_sat,
    drainsum)
    return drainsum

for compi in range (1):
# 1* calculate drainage of compartment


    delta_theta = Con(Compartment[compi].theta > Compartment[compi].theta_fc,
    calculate_delta_theta(Compartment[compi].theta, Compartment[compi].theta_sat, Compartment[compi].theta_fc, Compartment[compi].tau),
    0)
   
    delta_theta.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/adt')
   
    #drain_comp = delta_theta*Compartment[compi].thickness

# 2* check drainability

    excess = ZeroRaster
    pre_thick = ZeroRaster
    for i in range (compi):
        pre_thick = pre_thick + Compartment.thickness
    drainmax = delta_theta*pre_thick

    #drainmax.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/adm')

# 3* drain compartment

# 3.1: if drainsum <= drainmax

    drainsum_save = Con(drainsum<=drainmax,1,0)

    Compartment[compi].theta = Con(drainsum_save == 1,
    Compartment[compi].theta - delta_theta,
    Compartment[compi].theta)

    #Compartment[compi].theta.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/at')

    drainsum = Con(drainsum_save == 1,
    drainsum + drain_comp,
    drainsum)

    #drainsum.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/ads')

    excess = Con(drainsum_save == 1,
    CheckExcess(drainsum, excess,Compartment[compi].K_sat),
    excess)

    #excess.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/aex')
   
    drainsum = Con(drainsum_save == 1,
    CheckDrainsum(drainsum, Compartment[compi].K_sat),
    drainsum)

    #drainsum.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/ads1')

#3.2: if drainsum > drainmax
    delta_theta = Con(drainsum_save == 0,
    drainsum/pre_thick,
    delta_theta)

    #delta_theta.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/adt1')
   
    theta_x = Con(drainsum_save == 0,
    calculate_theta(delta_theta, Compartment[compi].theta_sat, Compartment[compi].theta_fc, Compartment[compi].tau),
    theta_x)

    #theta_x.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/atx')
   
#3.2.1: if theta_x <= theta_sat

    theta_x_save = Con(theta_x <= Compartment[compi].theta_sat,1, 0)
   
    Compartment[compi].theta = Con(drainsum_save == 0, Con(theta_x_save==1,
    Compartment[compi].theta + drainsum/Compartment[compi].thickness),
    Compartment[compi].theta)

    #Compartment[compi].theta.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/at1')
#3.2.1.1: if theta > theta_x

    theta_save = Con(Compartment[compi].theta > theta_x,1,0)
   
    drainsum = Con(drainsum_save == 0, Con(theta_x_save==1,Con(theta_save==1,
    ((Compartment[compi].theta - theta_x) + delta_theta)*Compartment[compi].thickness)),
    drainsum)

    #drainsum.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/ads2')
   
    excess = Con(drainsum_save == 0, Con(theta_x_save==1,Con(theta_save==1,
    CheckExcess(drainsum,excess,Compartment[compi].K_sat))),
    excess)

    #excess.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/aex1')
   
    drainsum = Con(drainsum_save == 0, Con(theta_x_save==1,Con(theta_save==1,
    CheckDrainsum(drainsum,Compartment[compi].K_sat))),
    drainsum)

    #drainsum.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/ads3')
   
    Compartment[compi].theta = Con(drainsum_save == 0, Con(theta_x_save==1,Con(theta_save==1,
    theta_x - delta_theta)),
    Compartment[compi].theta)

    #Compartment[compi].theta.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/at2')


drainsum.save('C:/Users/MattiasVO/Documents/ArcGIS/Default.gdb/at3')

[/PHP]