Spatially joins points with all zones of the NHD appending the COMID of the catchment that each point falls in

sjoinPointsNHD.py

import numpy as np
import geopandas as gpd
from geopandas.tools import sjoin

NHD_dir = "D:/NHDPlusV21"
inputs = np.load('%s/StreamCat_npy/zoneInputs.npy' % NHD_dir).item()
pt_file = "C:/Users/Rdebbout/Downloads/Stations_albers.shp"

pts = gpd.GeoDataFrame.from_file(pt_file)
l = "/NHDPlus"
for zone in inputs:
    hr = inputs[zone]
    print zone
    cats = gpd.GeoDataFrame.from_file("%s%s%s%s%s%sCatchment/Catchment.shp" % (NHD_dir,l,hr,l,zone,l))
    joined = sjoin(pts, cats, how="left", op="within")
    cols = joined.columns
    good = joined.ix[joined.FEATUREID.notnull()]
    good = good[cols +['FEATUREID']]
    good.columns = cols + ['COMID']