Porting Avenue code for Producing Building Shadows to ArcPy/Python for ArcGIS Desktop? [duplicate]

Question

Whuber provided an answer at Producing building shadows using ArcGIS Desktop? which required using Avenue code.

Any idea how to make it work in ArcGIS Desktop 10?

Answer 1

I've set this up as a ArcToolbox type of script, rather than field calculator as whuber did.

This is pretty much a straight port of whubers Avenue code.

EDIT: the script assumes the height is stored in a field in the featureclass attribute table, not a featureclass with 3D geometries (PolygonZ)

import arcpy,os,math
arcpy.env.overwriteOutput=True

buildings = arcpy.GetParameterAsText(0)
shadows = arcpy.GetParameterAsText(1)
heightfield=arcpy.GetParameterAsText(2) #Must be in the same units as the coordinate system!
shapefield=arcpy.Describe(buildings).shapeFieldName
try:azimuth=float(arcpy.GetParameterAsText(3))
except:azimuth=200 #default
try:altitude=float(arcpy.GetParameterAsText(4))
except:altitude=35 #default

#Output
result=arcpy.CreateFeatureclass_management(os.path.dirname(shadows),os.path.basename(shadows),'POLYGON')
inscur=arcpy.InsertCursor(shadows)

# Compute the shadow offsets.
spread = 1/math.tan(altitude) #outside loop as it only needs calculating once

for row in arcpy.SearchCursor(buildings):
    shape=row.getValue(shapefield)
    height=float(row.getValue(heightfield))

    # Compute the shadow offsets.
    x = -height * spread * math.sin(azimuth)
    y = -height * spread * math.cos(azimuth)

    # Clone the original shape.
    clone=arcpy.CopyFeatures_management(shape,arcpy.Geometry())[0]

    # Adjoin the wall shadows.
    for part in shape:
        for i,j in enumerate(range(1,part.count)):
            pnt0=part[i] #This will fail if the scripts comes across a polygon with
            pnt1=part[j] #inner ring/s, to handle this case you'll need to test
                         #that each point is not None.
            if pnt1 is None:break #EDIT: now it won't fail, but you still 
                                  #don't get nice shadows from the inner walls.
            pnt0offset=arcpy.Point(pnt0.X+x,pnt0.Y+y)
            pnt1offset=arcpy.Point(pnt1.X+x,pnt1.Y+y)
            arr=arcpy.Array([pnt0,pnt1,pnt1offset,pnt0offset,pnt0])
            clone=arcpy.Union_analysis([arcpy.Polygon(arr),clone],arcpy.Geometry())
            clone=arcpy.Dissolve_management(clone,arcpy.Geometry())[0]

    newrow=inscur.newRow()
    newrow.shape=clone
    inscur.insertRow(newrow)
    del newrow,clone,arr,pnt0,pnt0offset,pnt1,pnt1offset

del inscur

Answer 2

I made a few improvements to @Luke's version of the script, mainly for better performance. The GP methods really don't like being called in a tight loop; by removing an unnecessary GP call and moving the necessary one outside of the innermost loop I sped up performance by 10x at least.

This also fixes shadows not being created for inner rings, degrees are converted to radians, the original FID is stored as an attribute of the output feature for joining purposes, and progress is reported at 10% increments.

import arcpy
import os
import math

def message(msg, severity=0):
    # Adds a Message (in case this is run as a tool)
    # and also prints the message to the screen (standard output)
    #
    print msg

    # Split the message on \n first, so that if it's multiple lines,
    #  a GPMessage will be added for each line
    try:
        for string in msg.split('\n'):
            # Add appropriate geoprocessing message
            #
            if severity == 0:
                arcpy.AddMessage(string)
            elif severity == 1:
                arcpy.AddWarning(string)
            elif severity == 2:
                arcpy.AddError(string)
    except:
        pass

def main():
    arcpy.env.overwriteOutput=True

    # Read in parameters
    inputFC = arcpy.GetParameterAsText(0)
    outputFC = arcpy.GetParameterAsText(1)
    heightfield = arcpy.GetParameterAsText(2) #Must be in the same units as the coordinate system!
    azimuth = math.radians(float(arcpy.GetParameterAsText(3))) #Must be in degrees
    altitude = math.radians(float(arcpy.GetParameterAsText(4))) #Must be in degrees

    # Specify output field name for the original FID
    origfidfield = "ORIG_FID"

    # Get field names
    desc = arcpy.Describe(inputFC)
    shapefield = desc.shapeFieldName
    oidfield = desc.oidFieldName

    #Output
    message("Creating output feature class %s ..." % outputFC)
    arcpy.CreateFeatureclass_management(
        os.path.dirname(outputFC),
        os.path.basename(outputFC),
        'POLYGON', "", "", "",
        desc.spatialReference if not arcpy.env.outputCoordinateSystem else "")
    arcpy.AddField_management(outputFC, origfidfield, "LONG")
    inscur = arcpy.InsertCursor(outputFC)

    # Compute the shadow offsets.
    spread = 1/math.tan(altitude) #outside loop as it only needs calculating once

    count = int(arcpy.GetCount_management(inputFC).getOutput(0))
    message("Total features to process: %d" % count)

    searchFields = ",".join([heightfield, oidfield, shapefield])
    rows = arcpy.SearchCursor(inputFC, "", "", searchFields)

    interval = int(count/10.0) # Interval for reporting progress every 10% of rows

    # Create array for holding shadow polygon vertices
    arr = arcpy.Array()
    for r, row in enumerate(rows):
        pctComplete = int(round(float(r) / float(count) * 100.0))
        if r % interval == 0:
            message("%d%% complete" % pctComplete)
        oid = row.getValue(oidfield)
        shape = row.getValue(shapefield)
        height = float(row.getValue(heightfield))

        # Compute the shadow offsets.
        x = -height * spread * math.sin(azimuth)
        y = -height * spread * math.cos(azimuth)

        # Initialize a list of shadow polygons with the original shape as the first
        shadowpolys = [shape]

        # Compute the wall shadows and append them to the list
        for part in shape:
            for i,j in enumerate(range(1,part.count)):
                pnt0 = part[i]
                pnt1 = part[j]
                if pnt0 is None or pnt1 is None:
                    continue # skip null points so that inner wall shadows can also be computed

                # Compute the shadow offset points
                pnt0offset = arcpy.Point(pnt0.X+x,pnt0.Y+y)
                pnt1offset = arcpy.Point(pnt1.X+x,pnt1.Y+y)

                # Construct the shadow polygon and append it to the list
                [arr.add(pnt) for pnt in [pnt0,pnt1,pnt1offset,pnt0offset,pnt0]]
                shadowpolys.append(arcpy.Polygon(arr))
                arr.removeAll() # Clear the array so it can be reused

        # Dissolve the shadow polygons
        dissolved = arcpy.Dissolve_management(shadowpolys,arcpy.Geometry())[0]

        # Insert the dissolved feature into the output feature class
        newrow = inscur.newRow()
        newrow.setValue(origfidfield, oid) # Copy the original FID value to the new feature
        newrow.shape = dissolved
        inscur.insertRow(newrow)

        del row, newrow, shadowpolys, dissolved
    del inscur, rows

if __name__ == "__main__":
    main()

NOTE: Performance is still pretty bad: Roughly 1.25 hours for ~34k buildings on a fast machine, and something is leaking memory at an alarming rate -- roughly 1MB/sec -- it will max out the process's 2GB RAM limit fairly quickly at which point performance will grind to a near standstill. Using gc.collect() does not help, and it sticks around even after the script has completed, so I suspect it is the geometry objects not being released internally, as in this ESRI forum thread. To confirm this, one could adjust the script to use temporary feature classes on disk instead of geometry objects for the intermediate geometry storage.

I think that this algorithm could also benefit greatly from multiprocessing. Since the algorithm does not require features to be grouped, one could partition by rows instead of by geographic areas.

Answer 3

As alluded to in the comments of my previous answer (kept intact instead of edited, for comparison purposes), performance could be a lot better with additional optimizations (using the in_memory workspace instead of using geometry objects, move GP calls outside of loops where possible) as well as by utilizing multiprocessing.

Here is another version of the script that runs much faster, especially if you use all of your machine's processors. It uses the standard multiprocessing module in Python, using the Pool construct to queue up jobs that are partitioned by row ranges and processed in parallel by a set number of processes. It is configurable to use all, some, or 1 of your CPU processors, along with specifying the row partition lengths. See the comments in the "configuration" section for more details.

Basically if you tell it to use multiprocessing, it will partition the input feature class into ranges of OIDs, create multiple intermediate shapefiles, one per range, and finally merge them at the end. It also cleans up after itself by removing the intermediate shapefiles. There is some overhead to this, but the benefit of using all of your CPU power far outweighs the overhead in my testing (getting approximately 80% scaling efficiency with a dataset of 32k features).

Notes:

• Because it uses the in_memory workspace for the computationally intensive shadow calculations, take care to specify an appropriate partition size (procfeaturelimit) if you are sending in extremely large datasets. With the default settings the partition size will be the count of input rows divided by the number of cores, which is probably the most efficient configuration unless you start running out of memory.

• I have had zero luck using this as a script tool in ArcMap/ArcCatalog, even when set to run out of process. It appears that cursors are extremely slow when run out of process and called from a script tool, and subprocesses never seem to get created. So I am just running it from the PyScripter remote engine or at the Windows command line. If you can get it to work as a script tool, please let me know! Edit: It does work fine as a script tool if you configure it to not use multiprocessing (cores = 1, procfeaturelimit = 0) and run it in process.

---

import arcpy
import os
import math
import multiprocessing
import time

############################    Configuration:    ##############################
# Specify scratch workspace
scratchws = r"c:\temp\bldgshadowtest" # MUST be a folder, not a geodatabase!

# Specify output field name for the original FID
origfidfield = "ORIG_FID"

# Specify the number of processors (CPU cores) to use (0 to use all available)
cores = 0

# Specify per-process feature count limit, tune for optimal
# performance/memory utilization (0 for input row count divided by cores)
procfeaturelimit = 0

# TIP: Set 'cores' to 1 and 'procfeaturelimit' to 0 to avoid partitioning and
# multiprocessing completely
################################################################################

def message(msg, severity=0):
    print msg

    try:
        for string in msg.split('\n'):
            if severity == 0:
                arcpy.AddMessage(string)
            elif severity == 1:
                arcpy.AddWarning(string)
            elif severity == 2:
                arcpy.AddError(string)
    except:
        pass

def getOidRanges(inputFC, oidfield, count):
    oidranges = []
    if procfeaturelimit > 0:
        message("Partitioning row ID ranges ...")
        rows = arcpy.SearchCursor(inputFC, "", "", oidfield, "%s A" % oidfield)
        minoid = -1
        maxoid = -1
        for r, row in enumerate(rows):
            interval = r % procfeaturelimit
            if minoid < 0 and (interval == 0 or r == count - 1):
                minoid = row.getValue(oidfield)
            if maxoid < 0 and (interval == procfeaturelimit - 1 or r == count - 1):
                maxoid = row.getValue(oidfield)
            if minoid >= 0 and maxoid >= 0:
                oidranges.append([minoid, maxoid])
                minoid = -1
                maxoid = -1
        del row, rows
    return oidranges

def computeShadows(inputFC, outputFC, oidfield, shapefield, heightfield, azimuth, altitude, outputSR="", whereclause=""):
    # Set outputs to be overwritten just in case; each subprocess gets its own environment settings
    arcpy.env.overwriteOutput=True

    # Create in-memory feature class for holding the shadow polygons
    tempshadows = r"in_memory\tempshadows"
    arcpy.CreateFeatureclass_management(
        "in_memory",
        "tempshadows",
        "POLYGON", "", "", "",
        outputSR)
    arcpy.AddField_management(tempshadows, origfidfield, "LONG")

    # Open a cursor on the input feature class
    searchfields = ",".join([heightfield, oidfield, shapefield])
    rows = arcpy.SearchCursor(inputFC, whereclause, "", searchfields)

    # Open an insert cursor on the in-memory feature class
    tempinscur = arcpy.InsertCursor(tempshadows)

    # Create array for holding shadow polygon vertices
    arr = arcpy.Array()

    # Compute the shadow offsets.
    spread = 1/math.tan(altitude)

    # Read the input features
    for row in rows:
        oid = int(row.getValue(oidfield))
        shape = row.getValue(shapefield)
        height = float(row.getValue(heightfield))

        # Compute the shadow offsets.
        x = -height * spread * math.sin(azimuth)
        y = -height * spread * math.cos(azimuth)

        # Copy the original shape as a new feature
        tempnewrow = tempinscur.newRow()
        tempnewrow.setValue(origfidfield, oid) # Copy the original FID value to the new feature
        tempnewrow.shape = shape
        tempinscur.insertRow(tempnewrow)

        # Compute the wall shadow polygons and insert them into the in-memory feature class
        for part in shape:
            for i,j in enumerate(range(1,part.count)):
                pnt0 = part[i]
                pnt1 = part[j]
                if pnt0 is None or pnt1 is None:
                    continue # skip null points so that inner wall shadows can also be computed

                # Compute the shadow offset points
                pnt0offset = arcpy.Point(pnt0.X+x,pnt0.Y+y)
                pnt1offset = arcpy.Point(pnt1.X+x,pnt1.Y+y)

                # Construct the shadow polygon and insert it to the in-memory feature class
                [arr.add(pnt) for pnt in [pnt0,pnt1,pnt1offset,pnt0offset,pnt0]]
                tempnewrow.shape = arr
                tempnewrow.setValue(origfidfield, oid) # Copy the original FID value to the new feature
                tempinscur.insertRow(tempnewrow)
                arr.removeAll() # Clear the array so it can be reused

    # Clean up the insert cursor
    del tempnewrow, tempinscur

    # Dissolve the shadow polygons to the output feature class
    dissolved = arcpy.Dissolve_management(tempshadows, outputFC, origfidfield).getOutput(0)

    # Clean up the in-memory workspace
    arcpy.Delete_management("in_memory")

    return dissolved

if __name__ == "__main__":
    arcpy.env.overwriteOutput=True

    # Read in parameters
    inputFC = arcpy.GetParameterAsText(0)
    outputFC = arcpy.GetParameterAsText(1)
    heightfield = arcpy.GetParameterAsText(2) #Must be in the same units as the coordinate system!
    azimuth = math.radians(float(arcpy.GetParameterAsText(3))) #Must be in degrees
    altitude = math.radians(float(arcpy.GetParameterAsText(4))) #Must be in degrees

    # Get field names
    desc = arcpy.Describe(inputFC)
    shapefield = desc.shapeFieldName
    oidfield = desc.oidFieldName

    count = int(arcpy.GetCount_management(inputFC).getOutput(0))
    message("Total features to process: %d" % count)

    #Export output spatial reference to string so it can be pickled by multiprocessing
    if arcpy.env.outputCoordinateSystem:
        outputSR = arcpy.env.outputCoordinateSystem.exportToString()
    elif desc.spatialReference:
        outputSR = desc.spatialReference.exportToString()
    else:
        outputSR = ""

    # Configure partitioning
    if cores == 0:
        cores = multiprocessing.cpu_count()
    if cores > 1 and procfeaturelimit == 0:
        procfeaturelimit = int(math.ceil(float(count)/float(cores)))

     # Start timing
    start = time.clock()

    # Partition row ID ranges by the per-process feature limit
    oidranges = getOidRanges(inputFC, oidfield, count)

    if len(oidranges) > 0: # Use multiprocessing
        message("Computing shadow polygons; using multiprocessing (%d processes, %d jobs of %d maximum features each) ..." % (cores, len(oidranges), procfeaturelimit))

        # Create a Pool of subprocesses
        pool = multiprocessing.Pool(cores)
        jobs = []

        # Get the appropriately delmited field name for the OID field
        oidfielddelimited = arcpy.AddFieldDelimiters(inputFC, oidfield)

        # Ensure the scratch workspace folder exists
        if not os.path.exists(scratchws):
            os.mkdir(scratchws)

        for o, oidrange in enumerate(oidranges):
            # Build path to temporary output feature class (dissolved shadow polygons)
            # Named e.g. <scratchws>\dissolvedshadows0000.shp
            tmpoutput = os.path.join(scratchws, "%s%04d.shp" % ("dissolvedshadows", o))

            # Build a where clause for the given OID range
            whereclause = "%s >= %d AND %s <= %d" % (oidfielddelimited, oidrange[0], oidfielddelimited, oidrange[1])

            # Add the job to the multiprocessing pool asynchronously
            jobs.append(pool.apply_async(computeShadows, (inputFC, tmpoutput, oidfield, shapefield, heightfield, azimuth, altitude, outputSR, whereclause)))

        # Clean up worker pool; waits for all jobs to finish
        pool.close()
        pool.join()

         # Get the resulting outputs (paths to successfully computed dissolved shadow polygons)
        results = [job.get() for job in jobs]

        try:
            # Merge the temporary outputs
            message("Merging temporary outputs into output feature class %s ..." % outputFC)
            arcpy.Merge_management(results, outputFC)
        finally:
            # Clean up temporary data
            message("Deleting temporary data ...")
            for result in results:
                message("Deleting %s" % result)
                try:
                    arcpy.Delete_management(result)
                except:
                    pass
    else: # Use a single process
        message("Computing shadow polygons; using single processing ...")
        computeShadows(inputFC, outputFC, oidfield, shapefield, heightfield, azimuth, altitude, outputSR)

    # Stop timing and report duration
    end = time.clock()
    duration = end - start
    hours, remainder = divmod(duration, 3600)
    minutes, seconds = divmod(remainder, 60)
    message("Completed in %d:%d:%f" % (hours, minutes, seconds))

As for performance, using the same dataset as my previous answer that originally took 1.25 hours, this version completes in 8 minutes, 30 seconds using a single process and no partitioning, and 2 minutes 40 seconds using 4 processes and 4 partitions (jobs) of ~8k features each. RAM usage is also much more reasonable, using around 250MB for a single process, and around 150MB per process for multiprocessing.

For comparison purposes, here is a test dataset we can benchmark against. It's ~22k records from the City of San Francisco's building footprints dataset (~40k polygon parts totaling 1,076,060 vertices, including one duplicate closing vertex for each ring), and I kept only one field with the height in feet of the building (prsizeh). With 4 processes it took 5 minutes and 30 seconds to complete, and with a single process it took 17 minutes, 30 seconds. I think they are fairly complex geometries because they were converted from 3D models rather than digitized from aerial photographs. (6277 features have multiple parts. One has 36 parts, many of them slivers.)