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

closed as too broad by PolyGeo Jul 9 at 21:50

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.