¿Usar OGR y Shapely de manera más eficiente? [cerrado]

28

Estoy buscando algunas sugerencias sobre cómo hacer que mi código de Python sea más eficiente. Normalmente, la eficiencia no me importa, pero ahora estoy trabajando con un archivo de texto de ubicaciones de EE. UU. Con más de 1.5 millones de puntos. Con la configuración dada, toma aproximadamente 5 segundos ejecutar las operaciones en un punto; Necesito bajar esta figura hacia abajo.

Estoy usando tres paquetes GIS de Python diferentes para realizar algunas operaciones diferentes en los puntos y generar un nuevo archivo de texto delimitado.

  1. Utilizo OGR para leer un shapefile de límite de condado y obtener acceso a la geometría de límite.
  2. Shapely comprueba si un punto está dentro de alguno de estos condados.
  3. Si está dentro de uno, uso la biblioteca de archivos de formas de Python para extraer información de atributos del límite .dbf.
  4. Luego escribo información de ambas fuentes en un archivo de texto.

Sospecho que la ineficiencia radica en tener un bucle de 2 a 3 niveles ... no estoy muy seguro de qué hacer al respecto. Particularmente busco ayuda con alguien con experiencia en el uso de cualquiera de estos 3 paquetes, ya que es la primera vez que uso cualquiera de ellos.

import os, csv
from shapely.geometry import Point
from shapely.geometry import Polygon
from shapely.wkb import loads
from osgeo import ogr
import shapefile

pointFile = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New\NationalFile_20110404.txt"
shapeFolder = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New"
#historicBounds = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New\US_Counties_1860s_NAD"
historicBounds = "US_Counties_1860s_NAD"
writeFile = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New\NewNational_Gazet.txt"

#opens the point file, reads it as a delimited file, skips the first line
openPoints = open(pointFile, "r")
reader = csv.reader(openPoints, delimiter="|")
reader.next()

#opens the write file
openWriteFile = open(writeFile, "w")

#uses Python Shapefile Library to read attributes from .dbf
sf = shapefile.Reader("C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New\US_Counties_1860s_NAD.dbf")
records = sf.records()
print "Starting loop..."

#This will loop through the points in pointFile    
for row in reader:
    print row
    shpIndex = 0
    pointX = row[10]
    pointY = row[9]
    thePoint = Point(float(pointX), float(pointY))
    #This section uses OGR to read the geometry of the shapefile
    openShape = ogr.Open((str(historicBounds) + ".shp"))
    layers = openShape.GetLayerByName(historicBounds)
    #This section loops through the geometries, determines if the point is in a polygon
    for element in layers:
        geom = loads(element.GetGeometryRef().ExportToWkb())
        if geom.geom_type == "Polygon":
            if thePoint.within(geom) == True:
                print "!!!!!!!!!!!!! Found a Point Within Historic !!!!!!!!!!!!"
                print str(row[1]) + ", " + str(row[2]) + ", " + str(row[5]) + " County, " + str(row[3])
                print records[shpIndex]
                openWriteFile.write((str(row[0]) + "|" + str(row[1]) + "|" + str(row[2]) + "|" + str(row[5]) + "|" + str(row[3]) + "|" + str(row[9]) + "|" + str(row[10]) + "|" + str(records[shpIndex][3]) + "|" + str(records[shpIndex][9]) + "|\n"))
        if geom.geom_type == "MultiPolygon":
            for pol in geom:
                if thePoint.within(pol) == True:
                    print "!!!!!!!!!!!!!!!!! Found a Point Within MultiPolygon !!!!!!!!!!!!!!"
                    print str(row[1]) + ", " + str(row[2]) + ", " + str(row[5]) + " County, " + str(row[3])
                    print records[shpIndex]
                    openWriteFile.write((str(row[0]) + "|" + str(row[1]) + "|" + str(row[2]) + "|" + str(row[5]) + "|" + str(row[3]) + "|" + str(row[9]) + "|" + str(row[10]) + "|" + str(records[shpIndex][3]) + "|" + str(records[shpIndex][9]) + "|\n"))
        shpIndex = shpIndex + 1
    print "finished checking point"
    openShape = None
    layers = None


pointFile.close()
writeFile.close()
print "Done"
    
pregunta GrantD 20.05.2011 - 19:03

2 respuestas

21

El primer paso sería mover el shapefile abierto fuera del bucle de filas, está abriendo y cerrando el shapefile 1.5 millones de veces.

Para ser honesto, me gustaría meter todo el lote en PostGIS y hacerlo utilizando SQL en tablas indexadas.

    
respondido por el Ian Turton 20.05.2011 - 20:24
19

Un vistazo rápido a su código nos recuerda algunas optimizaciones:

  • Primero, compruebe cada punto contra el cuadro delimitador / sobre de los polígonos, para eliminar los valores atípicos obvios. Podría ir un paso más allá y contar el número de bboxes en que se encuentra un punto, si es exactamente uno, entonces no necesita ser probado contra la geometría más compleja (bueno, en realidad lo será si se encuentra en más más de uno, tendrá que probarse más. Podría hacer dos pasos para eliminar los casos simples de los casos complejos).

  • En lugar de recorrer cada punto y probarlo contra polígonos, recorre los polígonos y prueba cada punto. La carga / conversión de la geometría es lenta, por lo que debe hacerlo lo menos posible. Además, cree inicialmente una lista de puntos del CSV, para evitar tener que hacerlo varias veces por punto y luego descarte los resultados al final de esa iteración.

  • Indexe espacialmente sus puntos, lo que implica convertirlos en un shapefile, SpatialLite o algo así como un PostGIS / PostgreSQL de la base de datos. Esto tiene la ventaja de que herramientas como OGR podrán hacer la mayor parte del trabajo por usted.

  • No escriba la salida hasta el final: print () es una función costosa en el mejor de los casos. En su lugar, almacene los datos como una lista y escríbalos al final utilizando las funciones de decapado de Python o las funciones de volcado de listas.

respondido por el MerseyViking 20.05.2011 - 19:37

Lea otras preguntas en las etiquetas