Whuber proporcionó una respuesta en ¿Produciendo sombras de edificios utilizando ArcGIS Desktop? que requirió el uso del código de Avenue .
¿Alguna idea de cómo hacer que funcione en ArcGIS Desktop 10?
Whuber proporcionó una respuesta en ¿Produciendo sombras de edificios utilizando ArcGIS Desktop? que requirió el uso del código de Avenue .
¿Alguna idea de cómo hacer que funcione en ArcGIS Desktop 10?
He configurado esto como un tipo de script de ArcToolbox, en lugar de una calculadora de campo como lo hizo whuber.
Esto es más o menos un puerto directo de Whubers Avenue code.
EDITAR: el script asume que la altura se almacena en un campo en 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
Hice algunas mejoras en la versión de @ Luke del script, principalmente para un mejor rendimiento. A los métodos de GP realmente no les gusta que los llamen en un circuito cerrado; al eliminar una llamada GP innecesaria y mover la necesaria fuera del bucle más interno, aceleré el rendimiento al menos 10 veces.
Esto también corrige las sombras que no se crean para los anillos internos, los grados se convierten en radianes, el FID original se almacena como un atributo de la función de salida para propósitos de unión y el progreso se informa 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:elrendimientosiguesiendobastantemalo:Aproximadamente1.25horaspara~34kedificiosenunamáquinarápida,yalgoestáperdiendomemoriaaunavelocidadalarmante,aproximadamente1MB/seg.ellímitedememoriaRAMde2GBdelprocesoconbastanterapidez,momentoenelqueelrendimientosedetendrácasialmínimo.Elusodegc.collect()
noayuda,ysemantieneinclusodespuésdequeelscriptsehayacompletado,porloquesospechoquelosobjetosdegeometríanoselanzaninternamente,comoeneste
Creo que este algoritmo también podría beneficiarse enormemente de multiproceso . Dado que el algoritmo no requiere que las características se agrupen, se podría dividir por filas en lugar de por áreas geográficas.
Como se mencionó 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 (usando el espacio de trabajo in_memory
en lugar de usar objetos de geometría, mover llamadas de GP fuera de los bucles cuando sea posible), así como utilizando multiprocesamiento .
Aquí hay otra versión del script que se ejecuta mucho más rápido , especialmente si usa todos los procesadores de su máquina. Utiliza el módulo estándar multiprocessing
en Python, usando el Pool
para crear una cola de trabajos divididos por rangos de filas y Procesado en paralelo por un número determinado de procesos. Es configurable para usar todos, algunos o uno de sus procesadores de CPU, junto con la especificación de las longitudes de partición de la fila. Consulte los comentarios en la sección "configuración" para obtener más detalles.
Básicamente, si le dice que use el multiprocesamiento, dividirá la clase de entidad de entrada en rangos de OID, creará varios shapefiles intermedios, uno por rango, y finalmente los combinará al final. También se limpia después de sí mismo mediante la eliminación de los shapefiles intermedios. Esto conlleva cierta sobrecarga, pero el beneficio de usar toda la potencia de su CPU supera con creces la sobrecarga en mis pruebas (obteniendo una aproximadamente 80% de eficiencia de escala con un conjunto de datos de 32k características).
Notas :
Debido a que usa el espacio de trabajo in_memory
para los cálculos de sombra intensivos computacionalmente, tenga cuidado de especificar un tamaño de partición apropiado ( procfeaturelimit
) si está enviando conjuntos de datos extremadamente grandes. Con la configuración predeterminada, el tamaño de la partición será el número de filas de entrada dividido por el número de núcleos, que probablemente sea la configuración más eficiente a menos que empiece a quedarse sin memoria.
No he tenido suerte al usar esto como una herramienta de script en ArcMap / ArcCatalog , incluso cuando está configurado para quedarse sin proceso. Parece que los cursores son extremadamente lentos cuando se agotan los procesos y se llaman desde una herramienta de script, y parece que los subprocesos nunca se crean. Así que solo lo estoy ejecutando desde PyScripter motor remoto o en la línea de comando de Windows. Si puede hacerlo funcionar como una herramienta de script, ¡hágamelo saber! Editar: Funciona bien como herramienta de script si lo configuras para que no use el multiprocesamiento ( cores = 1
, procfeaturelimit = 0
) y lo ejecutes en el 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 , usando el mismo conjunto de datos que mi respuesta anterior que originalmente tomó 1.25 horas, esta versión se completa en 8 minutos, 30 segundos usando un solo 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 RAM también es mucho más razonable, ya que utiliza alrededor de 250 MB para un solo proceso y alrededor de 150 MB por proceso para multiprocesamiento.
A efectos de comparación , aquí hay un conjunto de datos de prueba que podemos comparar. Son ~ 22k registros de conjunto de datos de huellas de la Ciudad de San Francisco (~ 40k Partes de polígonos que totalizan 1,076,060 vértices, incluido un vértice de cierre duplicado para cada anillo), y mantuve solo un campo con la altura en pies del edificio ( prsizeh
). Con 4 procesos, tomó 5 minutos y 30 segundos para completarse, y con un solo proceso tomó 17 minutos, 30 segundos. Creo que son geometrías bastante complejas porque se convirtieron de modelos 3D en lugar de digitalizarse a partir de fotografías aéreas. (Las características del 6277 tienen varias partes. Una tiene 36 partes, muchas de ellas astilladas).
Lea otras preguntas en las etiquetas arcpy arcgis-10.0 parallel-processing