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]