Cómo encontrar las coordenadas de un CSV grande dentro de múltiples polígonos de shapefile

Question

Tengo un archivo CSV de unos 230 millones de filas, que contiene una columna de latitud, una columna de longitud, un sello de fecha y una columna de identificación. Se ve así (he añadido la cabecera en este ejemplo para ayudar a explicar, el archivo no tiene una cabecera):

ID,                               Date,                     Lon,     lat
46d4089a713082d85452a2af64571644, 2016-11-30T12:57:11.000Z, 53.4529, -2.287
46d4089a713082d85452a2af64571644, 2016-11-30T12:57:26.000Z, 53.4521, -2.2859
46d4089a713082d85452a2af64571644, 2016-11-30T12:57:59.000Z, 53.4522, -2.2878
46d4089a713082d85452a2af64571644, 2016-11-30T12:59:01.000Z, 53.4547, -2.284
a6af7b30dc3ffa0ee7ecea02a2981b7d, 2016-11-30T13:03:01.000Z, 53.4457, -2.2774
693316c8b95e9fb07207400414714180, 2016-11-30T21:18:16.913Z, 53.2887, -2.6687

También tengo dos shapefiles de polígonos ESRI que he creado. Lo que quiero es una lista de todos los ID que tienen al menos una coordenada en ambos shapefiles.

Como el archivo es tan grande y no podía cargarse en la memoria, mi enfoque fue dividir el archivo CSV en 117 CSV de 2 millones de filas. Luego planeé usar geopandas para leer los CSVs más pequeños como un GeodataFrame, y usar sjoin para encontrar todos los waypoints en el shapefile 1.

Entonces tomaría la columna ID de los waypoints dentro del shapefile como una lista, tomaría todos los waypoints en el csv grande con un ID en la lista, y haría lo mismo de nuevo encontrando los puntos en el shapefile 2.

He intentado utilizar la indexación espacial para acelerar el proceso, ya que es un archivo grande

Intenté probar esto usando el primer CSV de 2 millones de filas, usando este código:

import pandas as pd
import geopandas
from geopandas.tools import sjoin
from shapely.geometry import *

waypoints = pd.read_csv(r'largefile_1.csv',sep=',', names=['TripId','lat','lon'],usecols=[0,2,3])

waypoints['geometry'] = waypoints.apply(lambda x: Point((float(x.lon), float(x.lat))), axis=1)

point = geopandas.GeoDataFrame(waypoints, geometry='geometry')
polygon  = geopandas.GeoDataFrame.from_file(r'Shapefile1.shp')
point.crs = polygon.crs

spatial_index = point.sindex
possible_matches_index = list(spatial_index.intersection(polygon.bounds))
possible_matches = point.iloc[possible_matches_index]
precise_matches = possible_matches[possible_matches.intersects(polygon)]

¿Hay alguna forma más eficiente de hacerlo, en otro programa o con diferentes plugins?

Answer 1

Si tienes acceso a una base de datos postgres (si no es así yo lo descargaría e instalaría la extensión postgis)

Aquí hay un ejemplo de cómo uso python y psycopg2:

import psycopg2
import psycopg2.extras

class DBException(Exception):
    pass
class DB():
    def __init__(self, parent_widget):
        self.conn = None
        self.dbhost = 'localhost'
        self.dbname = 'dbname'# Name your db 
        self.dbuser = 'postgres' #If you dont have set it something else
        self.dbpass = 'postgres' #what ever your password is

    def _connect(self):
        """
        Connects to the database
        :return:
        """
        if self.conn is None:
            try:
                self.conn = psycopg2.connect(
                    host=self.dbhost,
                    database=self.dbname,
                    user=self.dbuser,
                    password=self.dbpass
                    )
            except psycopg2.OperationalError as e:
                raise DBException("Error connecting to database on '%s'. %s" %
                                  (self.dbhost, str(e)))

    def _close(self):
        """
        Closes the connection to the database
        :return:
        """
        if self.conn is not None:
            self.conn.close()
            self.conn = None

    def create_table(self, sql, tbl_name):
        """
        :param sql: text string with the sql statement
        :param tbl_name: text string
        :return:
        """
        self._connect()
        cur = self.conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
        cur.execute("DROP TABLE IF EXISTS " + str(tbl_name))
        cur.execute(sql)
        self.conn.commit()
        self._close()

    def insert_data(self, sql):
        """
        :param sql: text string with the sql statement
        :return:
        """
        self._connect()
        cur = self.conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
        #print(sql)
        cur.execute(sql)
        self.conn.commit()
        self._close()
DB = DB()
DB.create_table("create table large_table(id varchar(40) PRIMARY KEY, date_insterted timestamp, pos geometry(POINT, 4326))" )

with open('the_csv_file.csv', 'r') as f:
    all_lines = f.readlines() # U might have to do some buffer work here
    sql_row_count = 0
    sql = ""
    for row in all_lines:
        row = row.split(',')
        if sql_row_count > 10000:
            DB.insert_data(sql)
            sql = ""
            sql_row_count = 0

        sql += "Insert into large_table id, date_insterted, pos Values(" + row[0] + ", " + row[1]+", ST_PointFromText('POINT(" + row[2] + " " + row[3]) + ")',4326); "

Eso es un comienzo para meterlos todos en una base de datos postgres, luego puedes escribir preguntas en postgres como "select * from large_table where st_x(pos) > 45 and st_x <55"

Espero que eso ayude.