[GIS] Polygonize.py outputs into a single file or simple merge

cartographygdalogrtopography

Given an unique raster topographic file at start, I want to get vector polygons for elevations thresholds 1m, 200m, 500m, 1000m, 2000m into a single shp file. I first slice my raster into 5 (or more) rasters, then polygonize each, and then have to merge them.

How to get these outputs into a single shape file ? Either by sending them as layers directly into the same file, or by creating individual shp files to merge into one final shp.

Current workflow details: It seems I must do this because poligonize works on pixels with a same value (elevation) x.

1. Generate raster slices:
I used gdal_calc.py to generate kinds of vertical slices, with one .tif for each level :

gdal_calc.py -A crop.tif --outfile=level001.tif --calc="1*(A>0)"     --NoDataValue=0
gdal_calc.py -A crop.tif --outfile=level200.tif --calc="200*(A>200)" --NoDataValue=0
...
...
...
...

All conserved pixels in level001.tif have an elevation of 1. All conserved pixels in level200.tif have an elevation of 200m. Opened together in QGis it give slices like this:

enter image description here

2. Polygonize my rasters slices (files):

I run polygonize.py separatedly on each of these several raster files to get my vector layers:

gdal_polygonize.py level001.tif -f "ESRI Shapefile" levels001.shp levels001 elev
gdal_polygonize.py level200.tif -f "ESRI Shapefile" levels200.shp levels200 elev
...
...
...
...

3. Merge into one shp file.

How to merge these shp outputs ?

Alternatively, I'am also interested by smoother workflow (loop).

Best Answer

If the shapefiles to be merged are in the same directory:

ogr2ogr merge.shp levels001.shp
ogr2ogr -update -append merge.shp levels200.shp -nln merge
ogr2ogr -update -append merge.shp levels500.shp -nln merge
ogr2ogr -update -append merge.shp levels1000.shp -nln merge
ogr2ogr -update -append merge.shp levels2000.shp -nln merge

Related Solutions

[GIS] How to conditionnally assign a new value to pixels of a raster image

The best all round tool here is a raster calculator.

gdal_calc is a GDAL raster calculator implemented in Python here, with some examples here.

If you e.g. wants to keep values above +50:

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>50)" --NoDataValue=0

You can specify several files -A to -Z, where each of them get a corresponding variable in the expression.

In this case, the boolean expression (A>50) returns either 1 or 0, and is then multiplied with A. The result is then stored in the corresponding pixel in outputfile.

Not sure if you need a programmatic approach or a manual one, but QGis has a raster calculator. See related question with answer here.

[GIS] Overlapping Rasters as numpy arrays

Because you are limited to writing contiguous blocks of data to GDAL rasters, the best way around this is to create the output array first, then write the array to the output raster. Using this method, you can use numpy's boolean indexing to dictate what is written. Your last for loop would work as follows using this logic:

import numpy

no_data = -2147483647
output_array = numpy.full((newheight, newwidth), no_data, 'int32')

for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform()
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    xoffset = int((f_xmin - xmin) / f_pwidth)
    yoffset = int((f_ymax - ymax) / f_pheight)  # Assumes grids align
    band = raster.GetRasterBand(1)
    data = band.ReadAsArray()
    # Create a boolean array that allows data where they exist
    mask = (output_array == no_data) & (data != band.GetNoDataValue())  # Note: this assumes you prefer data in order based on the files list
    # Allocate data to output array
    output_array[yoffset:yoffset + rows, xoffset:xoffset + cols][mask] = data[mask]

# Write data to the output raster
outband.WriteArray(output_array)
outband.FlushCache()

However, you will need to fully load the raster in to memory prior to flushing it to the disk. If memory-management is an issue, you can either use numpy.memmap, or read the overlapping block from the output raster for each iteration, which would look like this:

for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform() #geotransform for the old files being iterated over
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    xoffset = int((f_xmin - xmin) / f_pwidth)
    yoffset = int((f_ymax - ymax) / f_pheight)
    band = raster.GetRasterBand(1)
    data = band.ReadAsArray(0, 0, cols, rows)
    data_update = outband.ReadAsArray(xoffset, yoffset, cols, rows)
    mask = (data_update == no_data) & (data != band.GetNoDataValue())
    data_update[mask] = data[mask]
    outband.WriteArray(data_update, xoffset, yoffset)
    outband.FlushCache()

One last general note of caution: your method assumes that raster grids and data types align.

Related Question