¿Utilizar OGR y Shapely de forma más eficiente?

Question

Estoy buscando algunas sugerencias sobre cómo hacer mi código python más eficiente. Normalmente la eficiencia no me importa, pero ahora estoy trabajando con un archivo de texto de ubicaciones de Estados Unidos con más de 1,5 millones de puntos. Con la configuración dada está tomando alrededor de 5 segundos para ejecutar las operaciones en un punto; Necesito bajar esta cifra.

Estoy utilizando tres paquetes SIG de Python diferentes para hacer algunas operaciones diferentes en los puntos y la salida de un nuevo archivo de texto delimitado.

1. Utilizo OGR para leer un shapefile de límites del condado y obtener acceso a la geometría de los límites.
2. Shapely comprueba si un punto está dentro de alguno de estos condados.
3. Si está dentro de uno, utilizo la Python Shapefile Library para extraer la información de atributos del .dbf de los límites.
4. A continuación, escribo la información de ambas fuentes en un archivo de texto.

Sospecho que la ineficacia radica en tener un bucle de 2-3 niveles... no sé muy bien qué hacer al respecto. En particular, estoy buscando la ayuda de alguien con experiencia en el uso de cualquiera de estos 3 paquetes, ya que es mi primera vez para utilizar 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"

Answer 1

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

Para ser honesto, yo metería todo en PostGIS y lo haría usando SQL en tablas indexadas.

Answer 2

Un rápido vistazo a su código me hace pensar en algunas optimizaciones:

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

• En lugar de hacer un bucle a través de cada punto y probarlo contra los polígonos, haga un bucle a través de los polígonos y pruebe cada punto. La carga/conversión de la geometría es lenta, por lo que querrá hacerlo lo menos posible. Además, cree una lista de puntos desde el CSV inicialmente, de nuevo para evitar tener que hacerlo varias veces por punto y luego descartar los resultados al final de esa iteración.

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

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