I'd like to suggest you a few clusterization algorithms to do this.
Of course for that you should use some python libraries for machine learning. One of the most popular is scikit learn: http://scikit-learn.org/stable/index.html
Here you have an overview of clustering methods. There are a lot of algorithms and each is different. But the result is pretty much the same. You can look for the algorithm that takes number of samples as an input (for your task it will be 643 to 653) and perform the clusterization. Samples are always near to each others, so this will be good result for you.
Other possibility is to take simple k-means alghorithm and calculate the number of clusters:
n = number_of_points/648
cls = sklearn.cluster.KMeans(n_clusters=n)
cls.fit(x) # x is a matrix with your points
centroids = cls.cluster_centers_ # coordinates of centers for each cluster
classification = cls.predict(x) # cluster number for each point
It should divide your set to equal clusters, so just add a new column to your shapefile with prediction and dissolve this layer with this attribute.
Of course first you have to extract only coordinates, perform clusterization on that, and join classified coordinates with the shapefile.