patternpythonMinor
Raster processing for climate change model
Viewed 0 times
climaterastermodelforprocessingchange
Problem
I am working on a project on how climate change may affect the development of spruce budworm larvae throughout each year (historical and projected) using different climate change models and emission scenarios.
The code that I have works, but would take a ridiculous amount of time to run: 38 years when extrapolated out. I'm using python 2.7.12 on a 64-bit windows machine with 24GB RAM and xeon processors. So, I don't think machine specs are an issue... it's how my script is written and organized.
I'm fairly new to Python, and have read through some articles on NumPy arrays and scipy packages, but do not have the knowledge/skills on how to implement them in this script to optimize the raster processing. Is numpy or scipy the way to go, and can those packages accomplish what I am trying to do?
Each time a new 4-h timestep raster is accessed, 27 conditional statements are executed, which probably is the reason the script runs so slow. Does anyone see a way to reduce the amount of conditional statements, but achieve the same goal?
Any help on how to improve the overall performance so that the script runs in a reasonable amount of time would be greatly appreciated!
```
import arcpy
from arcpy.sa import *
import math
arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
ClimateModels = ["CanESM2", "CSIRO-Mk3-6-0", "HadGEM2-ES"]
ConcentrationPathway = ["rcp26", "rcp45", "rcp85"]
#100 random numbers between 0.4 and 2.5, representing the stochastic growth rates of individual spruce budworm larvae ....need to be the same each time it loops through GCM+RCP combination
BudwormList = [2.1, 2.3, 0.6, 2.1, 1.7, 1.1, 1.6, 2.4, 2.1, 2.2, 0.9, 1.6, 0.6, 2.0, 1.7, 0.6, 0.7, 1.8, 1.3, 2.4,
2.1, 0.8, 2.2, 1.6, 0.5, 1.0, 1.6, 2.4, 2.5, 2.5, 0.8, 0.5, 2.5, 0.9, 1.2, 0.7, 1.8, 1.8, 0.4, 1.9,
0.8, 1.3, 0.6, 2.2, 0.7, 1.6, 1.0, 0.7, 0.8, 2.4, 1.7, 0.5, 1.1, 2.2, 0.9, 1.6, 1.8, 0.8, 2.3, 0.9,
0.8, 1.0, 0.5, 1.6, 2.5, 2.2, 2.0, 1.3, 1.6, 1.
The code that I have works, but would take a ridiculous amount of time to run: 38 years when extrapolated out. I'm using python 2.7.12 on a 64-bit windows machine with 24GB RAM and xeon processors. So, I don't think machine specs are an issue... it's how my script is written and organized.
I'm fairly new to Python, and have read through some articles on NumPy arrays and scipy packages, but do not have the knowledge/skills on how to implement them in this script to optimize the raster processing. Is numpy or scipy the way to go, and can those packages accomplish what I am trying to do?
Each time a new 4-h timestep raster is accessed, 27 conditional statements are executed, which probably is the reason the script runs so slow. Does anyone see a way to reduce the amount of conditional statements, but achieve the same goal?
Any help on how to improve the overall performance so that the script runs in a reasonable amount of time would be greatly appreciated!
```
import arcpy
from arcpy.sa import *
import math
arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
ClimateModels = ["CanESM2", "CSIRO-Mk3-6-0", "HadGEM2-ES"]
ConcentrationPathway = ["rcp26", "rcp45", "rcp85"]
#100 random numbers between 0.4 and 2.5, representing the stochastic growth rates of individual spruce budworm larvae ....need to be the same each time it loops through GCM+RCP combination
BudwormList = [2.1, 2.3, 0.6, 2.1, 1.7, 1.1, 1.6, 2.4, 2.1, 2.2, 0.9, 1.6, 0.6, 2.0, 1.7, 0.6, 0.7, 1.8, 1.3, 2.4,
2.1, 0.8, 2.2, 1.6, 0.5, 1.0, 1.6, 2.4, 2.5, 2.5, 0.8, 0.5, 2.5, 0.9, 1.2, 0.7, 1.8, 1.8, 0.4, 1.9,
0.8, 1.3, 0.6, 2.2, 0.7, 1.6, 1.0, 0.7, 0.8, 2.4, 1.7, 0.5, 1.1, 2.2, 0.9, 1.6, 1.8, 0.8, 2.3, 0.9,
0.8, 1.0, 0.5, 1.6, 2.5, 2.2, 2.0, 1.3, 1.6, 1.
Solution
First, fix the bug!
You're using
Now, I'll start my review with a couple of style guides. You can read more about them here.
-
You should put a blank line between each group of imports. More, absolute imports are recommended, as they are usually more readable and tend to be better behaved (or at least give better error messages) if the import system is incorrectly configured
-
use string formatting when you're concatenating strings: e.g:
With all of the above in mind, we'll have the following code (the formatting is a bit ugly because you have really long formulas but heh..):
```
import math
import sys
import arcpy
from arcpy.sa import Con, Raster
arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
climate_models = ["CanESM2", "CSIRO-Mk3-6-0", "HadGEM2-ES"]
concentration_pathway = ["rcp26", "rcp45", "rcp85"]
budworm_list = [2.1, 2.3, 0.6, 2.1, 1.7, 1.1, 1.6, 2.4, 2.1, 2.2, 0.9, 1.6, 0.6, 2.0, 1.7, 0.6, 0.7, 1.8, 1.3, 2.4,
2.1, 0.8, 2.2, 1.6, 0.5, 1.0, 1.6, 2.4, 2.5, 2.5, 0.8, 0.5, 2.5, 0.9, 1.2, 0.7, 1.8, 1.8, 0.4, 1.9,
0.8, 1.3, 0.6, 2.2, 0.7, 1.6, 1.0, 0.7, 0.8, 2.4, 1.7, 0.5, 1.1, 2.2, 0.9, 1.6, 1.8, 0.8, 2.3, 0.9,
0.8, 1.0, 0.5, 1.6, 2.5, 2.2, 2.0, 1.3, 1.6, 1.2, 2.4, 0.9, 1.1, 0.8, 2.1, 2.2, 1.5, 2.4, 2.5, 2.2,
1.3, 0.4, 1.5, 1.3, 1.7, 1.9, 1.8, 0.4, 0.4, 0.5, 1.5, 0.5, 0.6, 0.7, 1.6, 0.9, 0.9, 2.5, 1.4, 1.3]
time_step = ['4h', '8h', '12h', '16h', '20h', '24h']
try:
for gcm in climate_models:
for rcp in concentration_pathway:
budworm_count = 1
for budworm in budworm_list:
year = 1971
while year 2.5) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 2) & (dev_l2 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 3) & (dev_l3 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 4) & (dev_l4 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 5) & (dev_l5 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 6) & (dev_l6_male > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 7) & (dev_pupa_male > 1)), dev_stage + 1, dev_stage)
dev_stage_female = Con(((dev_stage_female == 1) & (dev_l6_female > 1)),
dev_stage_female + 1, dev_stage_female)
dev_stage_female = Con(((dev_stage_female == 2) & (dev_pupa_female > 1)),
dev_stage_female + 1, dev_stage_female)
l2_date = Con(((dev_l2o > 1) & (l2_date == 0)), l2_date + current_day, l2_date)
l3_date = Con(((dev_l2 > 1) & (l3_date == 0)), l3_date + current_day, l3_date)
l4_date = Con(((dev_l3 > 1) & (l4_date == 0)), l4_date + current_day, l4_date)
l5_date = Con(((dev_l4 > 1) & (l5_date == 0)), l5_date + current_day, l5_date)
l6_date = Con(((dev_l5 > 1) & (l6_date == 0)), l6_date + current_day, l6_date)
pupa_male_date = Con(((dev_l6_male > 1) & (pupa_male_date == 0)),
pupa_male_date + current_day, pupa_male_date)
adult_male_date = Con(((dev_pupa_male > 1) & (adult_male_date == 0)),
adult_male_date + current_day, adult_male_date)
pupa_female_date = Con(((dev_l6_female > 1) & (pupa_female_date == 0)),
pupa_female_date + current_day, pupa_female_date)
adult_female_date = Con(((dev_pupa_female > 1) & (adult_female_date == 0)),
adult_female_date + current_day, adult_female_date)
current_day = current_day + 1
l2_date.save(
"C:\\Users\\Robert\\Documents\\PCI
You're using
sys.exc_info()[-1].tb_lineno but you're not importing the sys module anywhere, so do that first !Now, I'll start my review with a couple of style guides. You can read more about them here.
- use snake_case convention for your variable names
-
imports should be grouped in the following order:- standard library imports
- related third party imports
- local application/library specific imports
You should put a blank line between each group of imports. More, absolute imports are recommended, as they are usually more readable and tend to be better behaved (or at least give better error messages) if the import system is incorrectly configured
- don't add redundant parentheses in conditional statements
- use the
print()function even if you're using Python 2.7.12
- add a space before and after each operator
- try to keep your lines no longer than 120 characters
-
use string formatting when you're concatenating strings: e.g:
print("Skiwi hates stacks {} times more than I do!".format(100))With all of the above in mind, we'll have the following code (the formatting is a bit ugly because you have really long formulas but heh..):
```
import math
import sys
import arcpy
from arcpy.sa import Con, Raster
arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
climate_models = ["CanESM2", "CSIRO-Mk3-6-0", "HadGEM2-ES"]
concentration_pathway = ["rcp26", "rcp45", "rcp85"]
budworm_list = [2.1, 2.3, 0.6, 2.1, 1.7, 1.1, 1.6, 2.4, 2.1, 2.2, 0.9, 1.6, 0.6, 2.0, 1.7, 0.6, 0.7, 1.8, 1.3, 2.4,
2.1, 0.8, 2.2, 1.6, 0.5, 1.0, 1.6, 2.4, 2.5, 2.5, 0.8, 0.5, 2.5, 0.9, 1.2, 0.7, 1.8, 1.8, 0.4, 1.9,
0.8, 1.3, 0.6, 2.2, 0.7, 1.6, 1.0, 0.7, 0.8, 2.4, 1.7, 0.5, 1.1, 2.2, 0.9, 1.6, 1.8, 0.8, 2.3, 0.9,
0.8, 1.0, 0.5, 1.6, 2.5, 2.2, 2.0, 1.3, 1.6, 1.2, 2.4, 0.9, 1.1, 0.8, 2.1, 2.2, 1.5, 2.4, 2.5, 2.2,
1.3, 0.4, 1.5, 1.3, 1.7, 1.9, 1.8, 0.4, 0.4, 0.5, 1.5, 0.5, 0.6, 0.7, 1.6, 0.9, 0.9, 2.5, 1.4, 1.3]
time_step = ['4h', '8h', '12h', '16h', '20h', '24h']
try:
for gcm in climate_models:
for rcp in concentration_pathway:
budworm_count = 1
for budworm in budworm_list:
year = 1971
while year 2.5) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 4.4) & (tempfile 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 2) & (dev_l2 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 3) & (dev_l3 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 4) & (dev_l4 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 5) & (dev_l5 > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 6) & (dev_l6_male > 1)), dev_stage + 1, dev_stage)
dev_stage = Con(((dev_stage == 7) & (dev_pupa_male > 1)), dev_stage + 1, dev_stage)
dev_stage_female = Con(((dev_stage_female == 1) & (dev_l6_female > 1)),
dev_stage_female + 1, dev_stage_female)
dev_stage_female = Con(((dev_stage_female == 2) & (dev_pupa_female > 1)),
dev_stage_female + 1, dev_stage_female)
l2_date = Con(((dev_l2o > 1) & (l2_date == 0)), l2_date + current_day, l2_date)
l3_date = Con(((dev_l2 > 1) & (l3_date == 0)), l3_date + current_day, l3_date)
l4_date = Con(((dev_l3 > 1) & (l4_date == 0)), l4_date + current_day, l4_date)
l5_date = Con(((dev_l4 > 1) & (l5_date == 0)), l5_date + current_day, l5_date)
l6_date = Con(((dev_l5 > 1) & (l6_date == 0)), l6_date + current_day, l6_date)
pupa_male_date = Con(((dev_l6_male > 1) & (pupa_male_date == 0)),
pupa_male_date + current_day, pupa_male_date)
adult_male_date = Con(((dev_pupa_male > 1) & (adult_male_date == 0)),
adult_male_date + current_day, adult_male_date)
pupa_female_date = Con(((dev_l6_female > 1) & (pupa_female_date == 0)),
pupa_female_date + current_day, pupa_female_date)
adult_female_date = Con(((dev_pupa_female > 1) & (adult_female_date == 0)),
adult_female_date + current_day, adult_female_date)
current_day = current_day + 1
l2_date.save(
"C:\\Users\\Robert\\Documents\\PCI
Code Snippets
print("Skiwi hates stacks {} times more than I do!".format(100))import math
import sys
import arcpy
from arcpy.sa import Con, Raster
arcpy.CheckOutExtension('spatial')
arcpy.env.overwriteOutput = True
climate_models = ["CanESM2", "CSIRO-Mk3-6-0", "HadGEM2-ES"]
concentration_pathway = ["rcp26", "rcp45", "rcp85"]
budworm_list = [2.1, 2.3, 0.6, 2.1, 1.7, 1.1, 1.6, 2.4, 2.1, 2.2, 0.9, 1.6, 0.6, 2.0, 1.7, 0.6, 0.7, 1.8, 1.3, 2.4,
2.1, 0.8, 2.2, 1.6, 0.5, 1.0, 1.6, 2.4, 2.5, 2.5, 0.8, 0.5, 2.5, 0.9, 1.2, 0.7, 1.8, 1.8, 0.4, 1.9,
0.8, 1.3, 0.6, 2.2, 0.7, 1.6, 1.0, 0.7, 0.8, 2.4, 1.7, 0.5, 1.1, 2.2, 0.9, 1.6, 1.8, 0.8, 2.3, 0.9,
0.8, 1.0, 0.5, 1.6, 2.5, 2.2, 2.0, 1.3, 1.6, 1.2, 2.4, 0.9, 1.1, 0.8, 2.1, 2.2, 1.5, 2.4, 2.5, 2.2,
1.3, 0.4, 1.5, 1.3, 1.7, 1.9, 1.8, 0.4, 0.4, 0.5, 1.5, 0.5, 0.6, 0.7, 1.6, 0.9, 0.9, 2.5, 1.4, 1.3]
time_step = ['4h', '8h', '12h', '16h', '20h', '24h']
try:
for gcm in climate_models:
for rcp in concentration_pathway:
budworm_count = 1
for budworm in budworm_list:
year = 1971
while year <= 2070:
if (year % 400 == 0) or ((year % 4 == 0) and (year % 100 != 0)):
days = 366
else:
days = 365
dev_stage = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_Stage.tif")
dev_stage_female = Raster(
"C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_Stage_female.tif")
dev_l2o = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L2o.tif")
dev_l2 = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L2.tif")
dev_l3 = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L3.tif")
dev_l4 = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L4.tif")
dev_l5 = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L5.tif")
dev_l6_male = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L6_male.tif")
dev_l6_female = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_L6_female.tif")
dev_pupa_male = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_pupa_male.tif")
dev_pupa_female = Raster(
"C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Dev_pupa_female.tif")
l2_date = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Jul_L2.tif")
l3_date = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Jul_L3.tif")
l4_date = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Jul_L4.tif")
l5_date = Raster("C:\\Users\\Robert\\Documents\\PCIC_ClimateData\\TEST\\Jul_L5.tif")
ldef budworm_list(lower_limit, upper_limit, decimals):
return [round(random.uniform(lower_limit, upper_limit), decimals) for _ in range(100)]budworm_list = generate_budworm_list(0.4, 2.5, 1)...
for budworm_count, budworm in enumerate(budworm_list, start=1):
...Context
StackExchange Code Review Q#155610, answer score: 3
Revisions (0)
No revisions yet.