n = 32
xdim = np.linspace(np.min(xSkyCoords), np.max(xSkyCoords), n + 1)
ydim = np.linspace(np.min(ySkyCoords), np.max(ySkyCoords), n + 1)
grid = np.zeros((n, n))

for i in range(len(xSkyCoords)):
    x_pos = 0
    y_pos = 0
    for j in range(len(xdim) - 1):
        if xSkyCoords[i] >= xdim[j] and xSkyCoords[i] <= xdim[j + 1]:
            x_pos = j
    for k in range(len(ydim) - 1):
        if ySkyCoords[i] >= ydim[k] and ySkyCoords[i] <= ydim[j + 1]:
            y_pos = k
    grid[(y_pos, x_pos)] += 1