Creating a shapefile from a numpy array

Question

I was following this tutorial: http://www.acgeospatial.co.uk/k-means-sentinel-2-python/ on how to generate an unsupervised classification for sentinel-2 images.

img_ds = gdal.Open('data/products/sentinel_bands.tif', gdal.GA_ReadOnly)
band = img_ds.GetRasterBand(2)
img = band.ReadAsArray()
X = img.reshape((-1,1))
k_means = cluster.KMeans(n_clusters=8)
k_means.fit(X)
X_cluster = k_means.labels_
X_cluster = X_cluster.reshape(img.shape)

This generates the classification as a numpy array

Now, I would like to transform that array to a shapefile. From my research I saw that you have to transform the array to a geodataframe with geopandas. I was trying to do that here:

X_cluster_dataframe = pd.DataFrame(data=X_cluster[1:, 1:],
                                   index=X_cluster[1:, 0],
                                   columns=X_cluster[0, 1:])
X _cluster_geodataframe = gp.GeoDataFrame(X_cluster_dataframe,
                                          geometry=gp.points_from_xy())

But I do not know how to specify the geodataframe's geometry.

I know that there are tools that can transform rasters to shapefiles, but the idea is to implement an efficient way to classificate big sentinel-2 imagenery and it would be better if I can generate the shapefile for the classification without fisrt having to create the raster.

Answer 1

You could save the classification array as a raster in memory (using gdal's MEM driver to create it) and then use gdal.Polygonize() function. Here is an example of how to use it.

If you definitely do not want to create a raster file (not even in memory) another approach would be using the rasterio.features.shapes() function. What this function does is:

Get shapes and values of connected regions in a dataset or array.

and what it returns is:

polygon, value – A pair of (polygon, value) for each feature found in the image. Polygons are GeoJSON-like dicts and the values are the associated value from the image, in the data type of the image.

You can get these polygons using your array as input instead of a raster file. Here is an answer doing exactly that.

Then, it is only a matter of saving the polygons to a shapefile.