Advice on taming a slow loop required for visualization of large GIS dataset
The nested loop compares each poly to each point, so the complexity is (Npoly * Npts).  

If the grid is regular, that is each square in a row has the same top and bottom coordinate, and each square in a column has the same left and right coordinate, then binary searches (bisect in std library) could be used to determine which grid row and column a point falls in.  This would be Npts * log Npoly.  I think this could be combined with .groupby() to eliminate any explicit loops.  I don't have time to post code now,  maybe tomorrow.

(It's tomorrow)

Better yet, use numpy.histogram2d.  On my laptop, it takes about 220ms to bin 2M points:

import numpy as np

min_lat, max_lat = 49.22, 49.30
num_lat_bins = 10

min_long, max_long = -123.23, -122.92
num_long_bins = 20

num_points = 2_000_000

# random coordinates for testing
lats  = np.random.uniform(min_lat,  max_lat,  num_points)
longs = np.random.uniform(min_long, max_long, num_points)


h, lat_edges, long_edges = np.histogram2d(lats, longs,
                                          bins=(10,20),
                                          range=[(min_lat,max_lat), (min_long, max_long)])  


Where, h is the histogram, or count in each grid square and lat_edges and 'long_edges are the grid lines.  They look like:

# h  (10 rows by 20 columns)
[[ 9982. 10126. 10124. ... 10218. 10140.  9844.]
 [ 9971. 10096. 10035. ... 10057.  9844.  9923.]
 [ 9940. 10036.  9982. ...  9966. 10034.  9872.]
 ... 
 [10025. 10011. 10027. ...  9870.  9942.  9985.]]

# lat_edges
[49.22 , 49.228, 49.236, 49.244, 49.252, 49.26 , 49.268, 49.276,
       49.284, 49.292, 49.3  ]

# long_edges
[-123.23  , -123.2145, -123.199 , -123.1835, -123.168 , -123.1525,
 -123.137 , -123.1215, -123.106 , -123.0905, -123.075 , -123.0595,
 -123.044 , -123.0285, -123.013 , -122.9975, -122.982 , -122.9665,
 -122.951 , -122.9355, -122.92  ]