[GIS] XYZ Clustering in Python and numpy arrays

clusteringnumpypython

I am looking for a method of clustering XYZ data in Python that works on numpy arrays. I know of a few sources, such as clusterpy and Pysal but have had little success with them as they seem to still be in the development stage and do not support numpy arrays.

EDIT:

My numpy arrays are of XY points with a Z column of a value at the point such as soil moisture or plant height. I want to spatially cluster the Z values into similar regions.

Best Answer

scikit-learn has an extensive clustering library with many different methods available. As a bonus scikit-learn is one of the best documented Python libraries I've seen. When working with 3d point clouds I've had a lot of success with DBSCAN for instance.

Alternately as @Fezter suggests above, scipy offers two different methods of clustring: k-means (vector-quantization) and hierarchical classification. While there isn't the number of methods available scipy is often a good place to start as the clustering methods don't tend to require as much setup.

Related Solutions

[GIS] How to utilize NumPy arrays to optimize big data geoprocessing

I think the crux of the question here is which tasks in your workflow are not really ArcGIS dependent? Obvious candidates include tabular and raster operations. If the data must start and end within a gdb or some other ESRI format, then you need to figure out how to minimize the cost of this reformat (i.e., minimize the number of round trips) or even justify it--simply might be too expensive to rationalize. Another tactic is to modify your workflow to use python-friendly data models earlier (for instance, how soon could you ditch vector polygons?).

To echo @gene, while numpy/scipy are really great, don't assume that these are the only approaches available. You can also use lists, sets, dictionaries as alternative structures (although @blah238's link is pretty clear about efficiency differentials), there are also generators, iterators, and all kinds of other great, fast, efficient tools for working these structures in python. Raymond Hettinger, one of the Python developers, has all kinds of great general Python content out there. This video is a nice example.

Also, to add onto @blah238's idea on multiplexed processing, if you're writing/executing within IPython (not just the "regular" python environment), you can use their "parallel" package for exploiting multiple cores. I'm no whiz with this stuff, but find it a bit higher-level/newbie-friendly than the multiprocessing stuff. Probably really just an issue of personal religion there, so take that with a grain of salt. There's a good overview at about it starting at 2:13:00 in this video. The whole video is great for IPython in general.

Raster Reclassify Using Python, GDAL, and Numpy

Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where:

# reclassification    
lista[np.where( lista < 200 )] = 1
lista[np.where((200 < lista) & (lista < 400)) ] = 2
lista[np.where((400 < lista) & (lista < 600)) ] = 3
lista[np.where((600 < lista) & (lista < 800)) ] = 4
lista[np.where( lista > 800 )] = 5

On an array of 6163 by 3537 pixels (41.6mb) the classification is done in 1.59 seconds, where it takes 12min 41s using the double for loop. In total just a speedup of 478x.

Bottomline, never use a double for loop using numpy

Best Answer

Related Solutions

[GIS] How to utilize NumPy arrays to optimize big data geoprocessing

Raster Reclassify Using Python, GDAL, and Numpy

Related Question