How to find coordinates from large CSV within multiple shapefile polygons

Question

I have a CSV file of around 230 million rows, which contains a latitude column, a longitude column, a Date stamp, and an ID column. It look like this (I added the header in this example to help explain, the file does not have a header):

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

I also have two ESRI polygon shapefiles I created. What I want is a list of all the IDs that have at least one coordinate in both of the shapefiles.

Since the file is so large and could not be loaded into memory, my approach was to break the CSV file into 117 CSVs of 2 million rows. I then planned on using geopandas to read the smaller CSVs as a GeodataFrame, and use sjoin to find all the waypoints in shapefile 1.

I would then take the ID column of the waypoints within the shapefile as a list, take all waypoints in the large csv with an ID in the list, and do the same again finding the points in shapefile 2.

I tried to use spatial indexing to speed up the process, as it is a large file

I tried to test this using the first CSV of 2 million rows, using this code:

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

Is there a more efficient way of doing this, in another program or with different plugins?

Answer 1

If you have access to a postgres database (if not I would download it and install the extension postgis)

Here is an example how I uses python and 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); "

That is a start to get them all into a postgres database, then you can write questions in postgres like "select * from large_table where st_x(pos) > 45 and st_x <55"

Hope that helps!