Work-around for cell-size output from rasterstats package -- can't get through idx == 37

zstats_30M.py

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from rasterstats import point_query
from datetime import datetime as dt

nlcd = 'L:/Priv/CORFiles/Geospatial_Library/Data/Project/StreamCat/LandscapeRasters/QAComplete/nlcd2006.tif'
pre = 'D:/NHDPlusV21/NHDPlusGL/NHDPlus04'
cat = gpd.read_file('%s/NHDPlusCatchment/Catchment.shp' % pre)
cat.to_crs({'init':'epsg:5070'}, inplace=True)
#################################################################################
def pt_stats(cat):
    final = {}
    for idx, row in cat.iterrows():
        print idx
        value = row.FEATUREID
        geom = row.geometry
        xs = np.arange(geom.bounds[0]+15,geom.bounds[2],30)
        ys = np.arange(geom.bounds[1]+15,geom.bounds[3],30)    
        pts = []
        for x in xs:
            for y in ys:
                pt = Point(x,y)
                if pt.within(geom):
                    pts.append(pt)
        b = gpd.GeoDataFrame({'join':[value]*len(pts)},geometry=pts,crs=cat.crs)
        out = point_query(b,nlcd)
        out = [k for k in out if k is not None] # remove NoData, could track here   
        g = np.array(out)
        f = np.rint(g).astype(int)
        stats = {}
        for stat in f:
            if not stat in stats.keys():
                stats[stat] = 1
            else:
                stats[stat] += 1
        final[value] = stats
    return final
#################################################################################

if __name__ == '__main__':
    start = dt.now()
    this = pt_stats(cat.copy())
    print dt.now() - start
    pd.DataFrame(this).T