Python and multiple calls to Hillshade function

857
2
03-24-2013 08:35 PM
brucesimonson
New Contributor III
Hi Folks,

New to these forums, and python.  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:  It looks like my post may have lost indents needed for Python.  If these indents don't display correctly, you can assume I got them right.  Sorry for the inconvenience.
Tags (2)
0 Kudos
2 Replies
curtvprice
MVP Esteemed Contributor
Bruce - what is your goal here? Have you looked into the Solar Radiation toolset?

Desktop Help 10.1: Understanding solar radiation analysis
0 Kudos
brucesimonson
New Contributor III
Bruce - what is your goal here? Have you looked into the Solar Radiation toolset?

Desktop Help 10.1: Understanding solar radiation analysis


Hi Curt,

Yep, we've looked at the Solar Analysis toolset.  It is quite interesting; in fact, we are considering enhancing it to incorporate re-radiation in the model.  The main (downside) issue with this toolset, is that it can only be run from one (point) location at a time, and it does take a (long) while to run the whole thing.  We are interested in sun exposure across the entire watersheds, so running solar analysis at each point isn't at all feasible.  We will run it for selected points, though.

The main goal is to see where ridgelines and other geographic features cast shadows across these watersheds.  We've checked, and the software tool is remarkably close, matching up very well with shadows in orthophotos, where the time of the image acquisition is well known.  I am actually quite surprised by the accuracies we are seeing.

It's really too bad that the Hillshade_3d function isn't designed to accept rasters of slope and aspect as parameters, and thereby avoid the unnecessary step to recalculate these with every call to the function.  Or, at least that's what I'm seeing.  I'd love to hear of a work-around.

Thanks for your question, though.  It is definitely relevant.

---

The other questions:  do you know how to generate pyramids from within my python program?  And also, the file locks aren't clearing; there must be a way to clear these from within python.

---

Thanks again.

Cheers, from Juneau, Alaska,
-Bruce
0 Kudos