Aggressive use of Hillshade function, using Python, for example

[brucesimonson]

Hi Folks,

This is cross posted to Python and Raster forums.

The point for this (Spatial Analyst, or 3D Analyst) is the Hillshade function appears to generate a slope and aspect for each call.  If so, this is a lot of overhead for my application, which needs to call this function approximately 50000 times.

Read on:

New to these forums.  I wish to write a program which aggressively calls the Hillshade function.  Specifically, for two important watersheds near Juneau, I want to call this function with "SHADOWS" enabled, for 5 minute intervals over an entire year.  For our latitude, this means that, for the period from sunrise to sunset, over an entire year, there will be over 50000 calls to this function.

I have a working program, but I am hoping it can be streamlined.  Also, I can't get the program to clear locks, and I would like to add a call to build pyramids, but I can't get either of these to work.

I have built an external file, which I access with the module csv, that contains solar positions, at 5 minute intervals, over an entire year.  These solar positions are the azimuth and altitude needed by the hilllshade function.

Here are four sample rows from this file:

HLEF_HS_20130101_0630_120x97_00x51,20130101,120.97,0.51
HLEF_HS_20130101_0635_121x78_01x33,20130101,121.78,1.33
HLEF_HS_20130101_0640_122x59_02x14,20130101,122.59,2.14
HLEF_HS_20130101_0645_123x41_02x95,20130101,123.41,2.95

field 1:  filename (to name the created hillshade), with date, time, azimuth, and altitude in filename
field 2:  date
field 3:  azimuth
field 4:  altitude

Here's my program, at least the essential part:

#### generate hillshades
import arcpy, csv, string, os, win32com.client, datetime as dt
from arcpy import env

arcpy.env.workspace = "E:\\solar\\solar.gdb"

# file with solar positions
f = open('E:\\solar_data\\solarpos.csv', 'rt')

# DEM
z = 'E:\\solar_data\\behlef120921'

try:
    reader = csv.reader(f)
 
    n=0
 
    t0=dt.datetime.now()
 
    for row in reader:
 
        n+=1
 
        filename = row[0]
        azimuth  = row[2]
        altitude = row[3]
        hs_pathfile = "E:\\solar\\solar.gdb\\"+filename
 
        t1=dt.datetime.now()
 
        print n, hs_pathfile, azimuth, altitude, (t1-t0).seconds

        # call hillshade function
        arcpy.HillShade_3d(z, hs_pathfile, azimuth, altitude, "SHADOWS", 1 )

        # arcpy.BuildPyramids_management(hs_pathfile)    # commented out, doesn't seem to work 
 
        # clear resulting hillshade from ArcMAP, so the TOC doesn't overload
        mxd = arcpy.mapping.MapDocument('CURRENT')
        df = arcpy.mapping.ListDataFrames(mxd)[0]
 
        for lyr in arcpy.mapping.ListLayers(mxd, '', df):
            arcpy.mapping.RemoveLayer(df, lyr)
 
finally:
    f.close()

####

Issues:

1)  The documentation on hillshade seems to indicate that slope and aspect are computed as part of the hillshade function.  This is an immense amount of overhead.  Seems like there should be a way to pass the slope and aspect to this function, rather than computing them for each call.  Yikes.  Each iteration runs about 100 seconds.  Any ideas here?

2)  Looking at my directory (E:\\solar\\solar.gdb), where the file geodatabase stores the results of these calls to hillshade, record locks are being held, even after the last section of code in the read loop, which clears each of these rasters from the TOC after it is generated.  How can I clear these locks from the Python program?

3)  I can't get the BuildPyramids function to work.  I'd like to do this, and perhaps compute statistics, as part of the program, right after each hillshade raster is created.

Lots going on here.  Any help would be greatly, greatly appreciated.

Cheers!
-Bruce

PS:  I am running from ArcMap's Python window (ArcGIS 10.0).  I'd run from a different environment (e.g. command line), but I can't get this to work at the moment.  The ArcMAP Python window is my working access point to Python.  Python 2.7.3.

PPS:  I will post this to the raster forum as well.

PPPS:  Looks like the Python code lost necessary indents in this post.  You can assume I got this right; sorry for the inconvenience.