¿Portar el código de Avenue para producir sombras de edificios a ArcPy/Python para ArcGIS Desktop?

Question

Whuber proporcionó una respuesta en ¿Producir sombras de edificios con ArcGIS Desktop? que requería el uso del código de la Avenida.

¿Alguna idea de cómo hacerlo funcionar en ArcGIS Desktop 10?

Answer 1

He configurado esto como un tipo de script ArcToolbox, en lugar de calculadora de campo como whuber hizo.

Esto es más o menos un puerto directo del código de Whubers Avenue.

EDIT: el script asume que la altura se almacena en un campo de la tabla de atributos de featureclass, no en una featureclass con geometrías 3D (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

He hecho algunas mejoras a la versión de @Luke del script, principalmente para mejorar el rendimiento. A los métodos de GP no les gusta ser llamados en un bucle cerrado; eliminando una llamada de GP innecesaria y moviendo la necesaria fuera del bucle más interno, aceleré el rendimiento por lo menos 10 veces.

Esto también corrige las sombras que no se crean para los anillos interiores, los grados se convierten en radianes, el FID original se almacena como un atributo de la característica de salida para fines de unión, y el progreso se reporta en incrementos del 10%.

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()

NOTA: El rendimiento sigue siendo bastante malo: Aproximadamente 1,25 horas para ~34k edificios en una máquina rápida, y algo está perdiendo memoria a un ritmo alarmante - aproximadamente 1MB/seg. Utilizando gc.collect() no ayuda, y se queda incluso después de que el script haya terminado, así que sospecho que son los objetos de geometría los que no se liberan internamente, como en esto Hilo del foro de ESRI . Para confirmar esto, se podría ajustar el script para utilizar clases de características temporales en el disco en lugar de objetos de geometría para el almacenamiento de la geometría intermedia.

Creo que este algoritmo también podría beneficiarse mucho de multiprocesamiento . Como el algoritmo no requiere que las características estén agrupadas, se podría hacer una partición por filas en lugar de por zonas geográficas.

Answer 3

Como se alude en los comentarios de mi respuesta anterior (mantenida intacta en lugar de editada, para fines de comparación), el rendimiento podría ser mucho mejor con optimizaciones adicionales (utilizando el in_memory en lugar de utilizar objetos geométricos, mover las llamadas a GP fuera de los bucles siempre que sea posible), así como utilizando multiprocesamiento .

Esta es otra versión del script que se ejecuta mucho más rápido especialmente si utiliza todos los procesadores de su máquina. Utiliza el estándar multiprocessing en Python, utilizando el módulo Pool para poner en cola trabajos que se dividen por rangos de filas y se procesan en paralelo por un número determinado de procesos. Se puede configurar para utilizar todos, algunos o 1 de los procesadores de la CPU, junto con la especificación de las longitudes de las particiones de las filas. Vea los comentarios en la sección "configuración" para más detalles.

Básicamente, si le dices que utilice el multiprocesamiento, dividirá la clase de características de entrada en rangos de OIDs, creará múltiples shapefiles intermedios, uno por rango, y finalmente los fusionará al final. También limpia después de sí mismo mediante la eliminación de los shapefiles intermedios. Esto conlleva cierta sobrecarga, pero el beneficio de utilizar toda la potencia de su CPU supera con creces la sobrecarga en mis pruebas (obteniendo aproximadamente 80% de eficiencia de escalado con un conjunto de datos de 32k características).

Notas:

• Debido a que utiliza el in_memory para los cálculos de sombra, que son muy intensos desde el punto de vista informático, tenga cuidado de especificar un tamaño de partición adecuado ( procfeaturelimit ) si envía conjuntos de datos extremadamente grandes. Con la configuración por defecto, el tamaño de la partición será el recuento de filas de entrada dividido por el número de núcleos, que es probablemente la configuración más eficiente a menos que empiece a quedarse sin memoria.

• He tenido cero suerte usando esto como una herramienta de script en ArcMap/ArcCatalog incluso cuando está configurado para que se agote el proceso. Parece que los cursores son extremadamente lentos cuando se ejecuta fuera del proceso y se llama desde una herramienta de script, y los subprocesos nunca parecen crearse. Así que sólo lo estoy ejecutando desde el PyScripter motor remoto o en la línea de comandos de Windows. Si consigues que funcione como herramienta de script, ¡házmelo saber! Editar: Funciona bien como herramienta de script si lo configuras para que no utilice el multiprocesamiento ( cores = 1 , procfeaturelimit = 0 ) y ejecutarlo en proceso.

---

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))

En cuanto al rendimiento Utilizando el mismo conjunto de datos que mi respuesta anterior, que originalmente tardó 1,25 horas, esta versión se completa en 8 minutos y 30 segundos utilizando un único proceso y sin partición, y 2 minutos 40 segundos utilizando 4 procesos y 4 particiones (trabajos) de ~8k características cada uno. El uso de la memoria RAM también es mucho más razonable, utilizando unos 250MB para un solo proceso, y unos 150MB por proceso para el multiproceso.

A efectos de comparación aquí hay un conjunto de datos de prueba con la que podemos compararnos. Son ~22k registros de la ciudad de San Francisco huellas de edificios (~40k partes de polígonos con un total de 1.076.060 vértices, incluyendo un vértice de cierre duplicado para cada anillo), y conservé sólo un campo con la altura en pies del edificio ( prsizeh ). Con 4 procesos, se ha necesitado 5 minutos y 30 segundos para completar, y con un solo proceso tardó 17 minutos, 30 segundos. Creo que son geometrías bastante complejas porque se convirtieron a partir de modelos 3D en lugar de digitalizarse a partir de fotografías aéreas. (6277 rasgos tienen múltiples partes. Uno tiene 36 partes, muchas de ellas astillas).