topic Re: Pythonic Code Structure in Python Questions

Pythonic Code Structure

PeterWilson — Sat, 11 Dec 2021 09:03:16 GMT

I'm busy rewriting my hydrological model python scripts into functions that are modular and found that my original python script runs twice as fast as my updated python script. The only reason I can find is that the original scripts outputs are saved in memory before being passed onto the next function. My new python script is writing and reading each output from disk. I might be incorrect in my understanding. If anyone can give me guidance on improving my new python script in both pythonic structure and performance.

Original Python Script

'''
Created on May 20, 2015


@author: PeterW
'''
# import system modules and site packages
import os
import time
import arcpy
import ArcHydroTools


# check out Spatial Analyst Extension
arcpy.CheckOutExtension("Spatial")


# set environment settings
arcpy.env.overwriteOutput = True


# set input and output arguments
raw = r"E:\Python\Temp\Model04\DEM\raw"
rasWs = r"E:\Python\Temp\Model04\Layers04"
outWs = r"E:\Python\Temp\Model04\Model04.gdb"


# Processing time
def hms_string(sec_elapsed):
    h = int(sec_elapsed / (60 * 60))
    m = int(sec_elapsed % (60 * 60) / 60)
    s = sec_elapsed % 60
    return "{}h:{:>02}m:{:>05.2f}s".format(h, m, s)


start_time1 = time.time()


# archydro variables
fill_sinks = os.path.join(rasWs, "fil")
flow_dir = os.path.join(rasWs, "fdr")
flow_acc = os.path.join(rasWs, "fac")
streams = os.path.join(rasWs, "str")
stream_seg = os.path.join(rasWs, "strlnk")
catchment_grid = os.path.join(rasWs, "cat")
catchment_poly = os.path.join(outWs, "Layers","Catchment")
drainage_line = os.path.join(outWs, "Layers", "DrainageLine")
adj_catch = os.path.join(outWs, "Layers", "AdjointCatchment")


try:
    # calculate the fill sinks
    arcpy.AddMessage("Processing Fill Sinks")
    ArcHydroTools.FillSinks(raw, fill_sinks)
    
    # calculate the flow direction
    arcpy.AddMessage("Processing Flow Direction")
    ArcHydroTools.FlowDirection(fill_sinks, flow_dir)
        
    # calculate the flow accumulation
    arcpy.AddMessage("Processing Flow Accumulation")
    ArcHydroTools.FlowAccumulation(flow_dir, flow_acc)
    
    # calculate the maximum flow accumulation
    arcpy.AddMessage("Processing Flow Accumulation Maximum")
    maxcellsResult = arcpy.GetRasterProperties_management(flow_acc, "MAXIMUM")
    maxcells = maxcellsResult.getOutput(0)
    print maxcells
    
    # calculate the stream threshold number of cells
    arcpy.AddMessage("Processing Stream Threshold")
    stream_threshold_numcells = (int(maxcells)*0.125/100)
    print stream_threshold_numcells
    
    # calculate the stream definition
    arcpy.AddMessage("Processing Stream Definition")
    ArcHydroTools.StreamDefinition(flow_acc, stream_threshold_numcells, streams)
    
    # calculate the stream segmentation
    arcpy.AddMessage("Processing Stream Segmentation")
    ArcHydroTools.StreamSegmentation(streams, flow_dir, stream_seg)
    
    # calculate the catchment grid delineation
    arcpy.AddMessage("Processing Catchment Grid Delineation")
    ArcHydroTools.CatchmentGridDelineation(flow_dir, stream_seg, catchment_grid)
    
    # calculate the catchment polygons from the catchment grid
    arcpy.AddMessage("Processing Catchment Polygons")
    ArcHydroTools.CatchmentPolyProcessing(catchment_grid, catchment_poly)
    
    # calculate the drainage lines from the stream segmentation grid
    arcpy.AddMessage("Processing DrainageLines")
    ArcHydroTools.DrainageLineProcessing(stream_seg, flow_dir, drainage_line)
    
    # calculate the adjoint catchment polygons
    arcpy.AddMessage("Processing Ajdoint Catchments")
    ArcHydroTools.AdjointCatchment(drainage_line, catchment_poly, adj_catch)
    
    arcpy.AddMessage("Completed Processing archydro Main Model")


except:
    print(arcpy.GetMessages(2))
    pass


# Determine the time take to process hydrological characteristics
end_time1 = time.time()
print ("It took {} to process hydrological characteristics".format(hms_string(end_time1 - start_time1)))


arcpy.CheckInExtension("Spatial")

New Python Script

'''
Created on Feb 24, 2016


@author: PeterW
'''
# import system modules and site packages
import time
from pathlib import Path
import arcpy
import ArcHydroTools


# check out extension
arcpy.CheckOutExtension("Spatial")


# set environment settings
arcpy.env.overwriteOutput = True


# set input and output arguments
dem = r"E:\Python\Temp\Model04\DEM\raw"
raster_workspace = r"E:\Python\Temp\Model04\Layers04"
fgdb = Path(r"E:\Python\Temp\Model04\Model04.gdb")


# Processing time
def hms_string(sec_elapsed):
    h = int(sec_elapsed / (60 * 60))
    m = int(sec_elapsed % (60 * 60) / 60)
    s = sec_elapsed % 60
    return "{}h:{:>02}m:{:>05.2f}s".format(h, m, s)


start_time1 = time.time()


# generate fill sinks grid
def  fill_sinks(dem, raster_workspace):
    fill_sinks = "{0}\\{1}".format(raster_workspace, "fil")
    arcpy.AddMessage("Processing Fill Sinks")
    ArcHydroTools.FillSinks(dem, fill_sinks)

fill_sinks(dem, raster_workspace)


# generate flow direction grid
def flow_direction(raster_workspace):
    fill_sinks = "{0}\\{1}".format(raster_workspace, "fil")
    flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
    arcpy.AddMessage("Processing Flow Direction")
    ArcHydroTools.FlowDirection(fill_sinks, flow_direction)
  
flow_direction(raster_workspace)


# generate flow accumulation grid
def flow_accumulation(raster_workspace):
    flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
    flow_accumulation = "{0}\\{1}".format(raster_workspace, "fac")
    arcpy.AddMessage("Processing Flow Accumulation")
    ArcHydroTools.FlowAccumulation(flow_direction, flow_accumulation)


flow_accumulation(raster_workspace)


# calculate stream threshold based on 0.5% of maximum flow accumulation
def stream_threshold(raster_workspace):
    flow_accumulation = "{0}\\{1}".format(raster_workspace, "fac")   
    arcpy.AddMessage("Processing Flow Stream Threshold")
    maxcellsResult = arcpy.GetRasterProperties_management(flow_accumulation, "MAXIMUM")
    maxcells = maxcellsResult.getOutput(0)
    arcpy.AddMessage("{} Maximum Cells".format(maxcells))
    stream_threshold_numcells = (int(maxcells)*0.5/100)
    arcpy.AddMessage("{} Stream Threshold".format(stream_threshold_numcells))
    return stream_threshold_numcells


numcells = stream_threshold(raster_workspace)


# generate the stream definition grid
def stream_definition(raster_workspace):
    flow_accumulation = "{0}\\{1}".format(raster_workspace, "fac")
    stream = "{0}\\{1}".format(raster_workspace, "str")
    arcpy.AddMessage("Processing Stream Definition")
    ArcHydroTools.StreamDefinition(flow_accumulation, numcells, stream)


stream_definition(raster_workspace)


# generate the stream segmentation grid
def stream_segmentation(raster_workspace):
    stream = "{0}\\{1}".format(raster_workspace, "str")
    flow_direction = "{0}\\{1}".format(raster_workspace, "fac")
    stream_link = "{0}\\{1}".format(raster_workspace, "strlnk")
    arcpy.AddMessage("Processing Stream Segmentation")
    ArcHydroTools.StreamSegmentation(stream, flow_direction, stream_link)


stream_segmentation(raster_workspace)


# calculate the catchment grid delineation
def catchment_grid(raster_workspace):
    flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
    stream_link = "{0}\\{1}".format(raster_workspace, "strlnk")
    catchment_grid = "{0}\\{1}".format(raster_workspace, "cat")
    arcpy.AddMessage("Processing Catchment Grid Delineation")
    ArcHydroTools.CatchmentGridDelineation(flow_direction, stream_link, catchment_grid)


catchment_grid(raster_workspace)


# calculate the catchment polygons from the catchment grid
def catchment_polygon(raster_workspace, fgdb):
    catchment_grid = "{0}\\{1}".format(raster_workspace, "cat")
    catchment_polygon = "{0}\\{1}".format(Path(fgdb, "Layers"), "Catchment")
    arcpy.AddMessage("Processing Catchment Polygon")
    ArcHydroTools.CatchmentPolyProcessing(catchment_grid, catchment_polygon)


catchment_polygon(raster_workspace, fgdb)


# calculate the drainage lines from the stream segmentation grid
def drainage_line(raster_workspace, fgdb):
    stream_link = "{0}\\{1}".format(raster_workspace, "strlnk")
    flow_direction = "{0}\\{1}".format(raster_workspace, "fdr")
    drainage_line = "{0}\\{1}".format(Path(fgdb, "Layers"), "DrainageLine")
    arcpy.AddMessage("Processing DrainageLine")
    ArcHydroTools.DrainageLineProcessing(stream_link, flow_direction, drainage_line)


drainage_line(raster_workspace, fgdb)


# calculate the adjoint catchment polygons
def adjoint_catchment(fgdb):
    drainage_line = "{0}\\{1}".format(Path(fgdb, "Layers"), "DrainageLine")
    catchment_polygon = "{0}\\{1}".format(Path(fgdb, "Layers"), "Catchment")
    adjoint_catchment = "{0}\\{1}".format(Path(fgdb, "Layers"), "AdjointCatchment")
    arcpy.AddMessage("Processing Adjoint Catchment")
    ArcHydroTools.AdjointCatchment(drainage_line, catchment_polygon, adjoint_catchment)


adjoint_catchment(fgdb)


arcpy.CheckInExtension("Spatial") 


# Determine the time take to process hydrological characteristics
end_time1 = time.time()
print ("It took {} to process hydrological characteristics".format(hms_string(end_time1 - start_time1)))

Re: Pythonic Code Structure

DanPatterson_Retired — Sat, 11 Dec 2021 09:03:19 GMT

try timing the individual components to see where your bottleneck is. This decorator needs your functions to return something... even if it is nothing. Your functions don't return anything

from functools import wraps
#------- decorator -----------------------
def delta_time(func):
    """simple timing decorator function"""
    @wraps(func)
    def wrapper(*args, **kwargs):
        t0 = time.time()                # start time
        result = func(*args, **kwargs)  # ... run the function ...
        t1 = time.time()                # end time
        print("Results for... {}".format(func.__name__))
        print("  time taken ...{:12.9f} sec.".format(t1-t0))
        #print("\n {}".format(result))  # uncomment to print within wrapper
        return result                  # return the result of the function
    return wrapper

useage... just put the @delta_time line above functions that you want to time

@delta_time  
def alpha(input=0):  
    """useless function 1 """  
    return input

Re: Pythonic Code Structure

PeterWilson — Thu, 25 Feb 2016 13:32:27 GMT

Hi Dan

Thanks for the following, I'll have a look at the following this evening and get back to you if I get stuck, if that fine with you.

Re: Pythonic Code Structure

LukeWebb — Thu, 25 Feb 2016 17:07:49 GMT

Totally unrelated to your problem, but I noticed you are using string processing to build file paths. There is a built in python module that handles this for you, and works across operating systems etc if this is ever a concern.

In simple terms, it joins parameters, adding a slash in between.

fill_sinks = "{0}\\{1}".format(raster_workspace, "fil")

Becomes:

fill_sinks = os.path.join(raster_workspace, "fil")

You can also do as many parameters as you like to build longer paths.

fill_sinks = os.path.join(raster_workspace, "fil", "output", "test.shp")

Returns: 'E:\Python\Temp\Model04\Layers04\fil\output\test.shp'

Re: Pythonic Code Structure

PeterWilson — Thu, 25 Feb 2016 18:06:27 GMT

Hi Luke

If you look at my original Python Script I used os.path.join and decided to move to PathLib Path. The unfortunate thing is that all ArcPy scripts require paths to be string and not objects hence using strings to create the paths from the Path objects.

Re: Pythonic Code Structure

DanPatterson_Retired — Thu, 25 Feb 2016 18:33:06 GMT

Peter, I suspect most are still using python 2.7.x. For those that are interested, pathlib became available in python 3.4 more documentation is given here

11.1. pathlib — Object-oriented filesystem paths — Python 3.5.1 documentation

Re: Pythonic Code Structure

JasonTipton — Sat, 11 Dec 2021 09:03:22 GMT

Yeah.... I'm still on 2.7.x. What's the problem with os.path.join at 3.4? Not having 3.7 installed, but just looking through the documentation and playing around with a 3.x shell online, it seems that os.path.join returns a string object as it always has. It looks like the new pathlib is a way to store paths as objects instead of as strings so instead of having to do os.path.xyz(), you now may be able to do path.xyz.

Where I'm confused is how this affects the script. At the end of the day, Peter Wilson is just formatting the path object back to a string. I don't see why os.path.join wouldn't have worked here (not that it's better than using the pathlib). Also it looks path has its own ways to concatenate paths.

I went ahead and installed pathlib into 2.7. It looks like you can do this:

# constructor
p = Path('a', 'b', 'c')
print(p)
# append paths to end with the "/"
p = p / 'd'
print(p)
print(str(p))

and it will print

WindowsPath('a/b/c)
WindowsPath('a/b/c/d')
a\b\c

Either way

os.path.join
str(Path)      # and 
"{}".format(Path)

should all get you to the same end goal: A string representing a path.

Re: Pythonic Code Structure

DanPatterson_Retired — Thu, 25 Feb 2016 21:25:13 GMT

Jason

It is an alternate... that is all, the initial question was why Peter was doing it... the answer is because he is using python 3.4 and for his environment, he must prefer it.

I haven't completed my tests to see how it handles escape characters should one forget to use raw notation (there is a current thread on paths and path constructs and Darren posted a poll). One can also use glob as an alternate to os.

I think pathlib may be an attempt to further unify file handling under multiple platforms and environments. Further testing might reveal advantages that I haven't come across yet.

ADDENDUM

Filenames and file paths in Python

for addition quasi-related materials

Re: Pythonic Code Structure

JasonTipton — Thu, 25 Feb 2016 21:37:44 GMT

Group all of your methods at the top under imports instead of co-mingling functions and statements.
Define your variables once instead of multiple times. (flow_accumulation, flow_direction, etc) assuming they are always the same in every method.
- Line 98 (Script 2): It looks like you assigned flow accumulation "fac" to the flow_direction variable
Don't pass RasterWorkspace into every function. Declare it once and use it inside each method.
You don't really need methods.
- After you do the previous, the way it is currently written, you really just need print statements in between each statement letting the user know what is going on. Cutting back on recreating the same variable over and over will cut a lot of code out.

These won't improve performance*, but will make it more readable, less code to manage, and less code to mess up.

*It may insignificantly improve by nanoseconds by not having to re-create the same variable over and over again!

Re: Pythonic Code Structure

JasonTipton — Sat, 11 Dec 2021 09:03:25 GMT

I misread his comment. I thought that he had moved from os.path.join to Path because he was having issues with os.path.join() and arcpy. I now see that he moved *in favor* of Path but that arcpy isn't expecting a Path object but a string.

From a python perspective the Path object is an improvement because you are now dealing with a smarter object than a dumb string. The string doesn't know it is a path so you have to call some function to get the file name or directory or to append a subdirectory to the string. From looking at the methods, it looks like you can even create a directory or check if a path exists or not, all from the object itself without having to call another function.

As for escape characters, it looks like it doesn't handle them unless you handle them: same as a string or os.path.join(). At the end of the day, I think you're still better off letting python add your slashes and just pass in your directory names in a list and let python make your path for you, whether that be os.path.join(a,b,c,d,e,...) or Path(a,b,c) / d / e / ...

# WRONG -- Don't use \t or \a
p = Path("dog\turtle\ant")
p
>> WindowsPath('dog\turtle\x07nt')
print(p)
>> dog    urtlent

# RIGHT
p = Path("dog", "turtle", "ant")
p
>> WindowsPath('dog/turtle/ant')
print(p)
>> dog\turtle\ant

# If you need to append, you are using the / operator, not the literal character '/'
p / 'cat' / 'mouse
>> WindowsPath('dog/turtle/ant/cat/mouse')

Re: Pythonic Code Structure

DanPatterson_Retired — Thu, 25 Feb 2016 22:41:45 GMT

Jason...

Good, this will be useful to know for people working in ArcGIS Pro which is using python 3.4.x. ArcMap 10.4 is sticking with the python 2.7.x version.

So as people ... (transition from/stick with )... (ArcMap/ArcGIS Pro) things with become ..(clearer/more confusing).

Re: Pythonic Code Structure

PeterWilson — Fri, 26 Feb 2016 20:29:20 GMT

Hi Jason

Thanks for your feedback, its truly appreciated as I'm still a novice when it comes to programming. With regards to your suggestions. Would you mind giving me an example to explain your first suggestion:

Group all of your methods at the top under imports instead of co-mingling functions and statements.

I only realized once you pointed it out that I was defining my variables multiple times, thanks for pointing it out.

Lastly could you give me some advice how I could create a function that would build my input and output arguments based on the raster workspace and fgdb.

Re: Pythonic Code Structure

DarylVan_Dyke — Fri, 26 Feb 2016 20:46:00 GMT

Hi Peter-

Your observations about read/write speeds are really important for effective geoprocessing. Since you running in a stand-alone script, you can use the "IN_MEMORY" workspace to speed up many transactions. I'll let you read up on it, but essentially you create a virtual GDB called "IN_MEMORY" and set your environment workspace to reference it. For vector operations you can see huge improvements - raster as well, if you have the memory to support it.

More generally, I'd look at the read-write configurations on your workstation. Reading and writing to a server will be slow, reading and writing to a local SSD will be fast. You will also see improvements if you can separate out read and write media - e.g. read from one SSD, write to another. If you have lots of RAM, consider defining a RAM disk at the OS level - this will be virtual hard drive that occupies RAM, and will be Very, Very fast. But not as much room (typically ~4 GB for the free RamDisk tools).

As mentioned above, there are system profiling tools in Python to diagnose the relative cost of different implementations. You can see dramatic changes in Python processing time by changing from C-based/C-optimized intrinsic functions, to true-Python functions. Again, you can read up on that (especially with respect to iterators and memory models for hash tables).

Have fun!

Daryl

Re: Pythonic Code Structure

BrianGelder — Fri, 26 Feb 2016 21:37:03 GMT

I'll second Daryl's comments and add that 'IN_MEMORY' can also vastly speed up table operations such as joins, especially when both datasets are 'IN_MEMORY'. The improvements I have documented depend on the tool, but 2-20x improvements are quite common. However, for table manipulation it's still not nearly as fast as storing data in a true SQL database and then doing your SQL operations there. File servers are notoriously slow but I haven't tried splitting my reads/writes between SSDs or using RAMdisk. Thanks for the ideas!

Brian

Re: Pythonic Code Structure

PeterWilson — Sat, 27 Feb 2016 18:05:30 GMT

Hi Jason

After fixing the line 98 incorrect allocation of fac (Flow Accumulation) to Flow Direction, my script completes within the same amount of time as my previous script. I'm going to look into developing a function to allocate the variables automatically as this is standardised. The only input required from the user is the:

Input DEM
Raster Folder
File Geodatabase

Re: Pythonic Code Structure

JasonTipton — Sat, 11 Dec 2021 09:03:27 GMT

# imports at top
import os
import arcpy
# etc. . .

# define methods here
def fill_sinks():
    # fill sink code

def flow_direction():
    # flow direction code

def flow_accumulation():
    # flow accumulation code

# etc...

# write your statements inline here, and better yet test for "__main__" first
# to only execute this code if you are executing this file

if __name__ == '__main__'
    fill_sinks()
    flow_direction()
    flow_accumulation()
     # etc . . .

Re: Pythonic Code Structure

JasonTipton — Sat, 11 Dec 2021 09:03:30 GMT

Since this isn't a class, only a script, you really don't need to put this in a method. How you did it in your original script would work just fine:

# archydro variables 
fill_sinks = os.path.join(rasWs, "fil") 
flow_dir = os.path.join(rasWs, "fdr") 
flow_acc = os.path.join(rasWs, "fac") 
streams = os.path.join(rasWs, "str") 
stream_seg = os.path.join(rasWs, "strlnk") 
catchment_grid = os.path.join(rasWs, "cat") 
catchment_poly = os.path.join(outWs, "Layers","Catchment") 
drainage_line = os.path.join(outWs, "Layers", "DrainageLine") 
adj_catch = os.path.join(outWs, "Layers", "AdjointCatchment")

If you prefer Path, just replace os.path.join(x,y,z) with str(Path(x,y,z)).

Honestly, for what you are doing, I think you might have over-engineered your solution on the 2nd script. The 1st one didn't look too bad. For performance, make sure everything is working in memory like other people have suggested. Only output it to disk when you are through processing it.