3

My goal is convert a 2D numpy array to shapefile. The matrix results from a classification, so the numbers have a meaning.

I solved my problem, but the solution is somehow convoluted and complicated.

Does anyone has a better implementation?

The actual code is too long lines (you cand find it here implemented as a function array -> shapefile) for a relatively simple action.

One of the big problem is gdal_polygonize giving separate polygons for disjoint regions.

gdal_calc.py* seems a pretty good solution, but it isn't in the osgeo module.

Here an example of dummy data from (Raster Layers — Python GDAL/OGR Cookbook 1.0 documentation)*

plot_matrix = np.array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 1, 1, 1, 1, 0, 2, 2, 2, 2, 0, 1, 1, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 0, 0, 0, 0, 0, 2, 0, 2, 0, 1, 0, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 0, 1, 1, 0, 0, 2, 0, 2, 0, 1, 1, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 0, 0, 1, 0, 0, 2, 0, 2, 0, 1, 0, 1, 0, 2, 0, 0, 0],
                    [ 0, 1, 1, 1, 1, 0, 2, 2, 2, 2, 0, 1, 0, 1, 0, 2, 2, 2, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                    [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

Show the data

The data represent 2 disjoint region, with 2 different values of the classification (1 and 2) plus a background of "No Data" values, represented by 0s.

My actual solutions are :

Workflow #1 (Write some files around)

  1. write the 2D numpy array to GeoTiff

GeoTiff

  1. use (gdal_polygonize.py)* to convert to "vector polygons for all connected regions of pixels in the raster sharing a common pixel value" (see documentation). This produce 5 different polygons, for the 5 regions in my data.
  2. open the produced shapefile and aggregate the polygon sharing the same value in multi polygon collection, in my case there will be 2 MultiPolygon (0 represent no data in the original data)

Workflow #2 (everything in memory)

  1. The same as above, but all driver in 1 is 'MEM' , so no temporary .tiff is created.
  2. Also done in memory (ogr driver 'Memory' )
  3. Same as above.

The result is a shapefile with 2 MultiPolygon, each with an unique value of the Field "class"

enter image description here

I need at least 10 reputation to post more than 2 links, so the references are all here:

  • Raster Layers — Python GDAL/OGR Cookbook 1.0 documentation : pcjericks.github.io/py
  • gdalogr-cookbook/raster_layers.html#create-raster-from-array
  • gdal_polygonize.py: www.gdal.org/gdal_polygonize.html
  • gdal_calc.py : www.gdal.org/gdal_calc.html