i did a classification with a min/max threshold for stormdamaged areas. For that i calculated two indices (NSDVI) and (SAVI). I intersected two numpy arrays from those results. I recieve a binary GeoTIFF file. For my end result i want to have a shapefile (GeoJSON or ESRI). This is my code so far:

import numpy, sys
from osgeo import gdal
ds = gdal.Open('H:/Threshold/nsdvi.tif')
ds1 = gdal.Open('H:/Threshold/savi.tif')
rows = ds.RasterYSize
cols = ds.RasterXSize
myarray = numpy.array(ds.GetRasterBand(1).ReadAsArray())
myarray1 = numpy.array(ds1.GetRasterBand(1).ReadAsArray())
selection1 = numpy.logical_and(myarray1 < 0.4032,myarray1 > 0.0803)
selection = numpy.logical_and(myarray < 0.0005 ,myarray >-0.0255)
new2 = numpy.logical_and(selection==1,selection1==1)

driver = ds.GetDriver()
outDs = driver.Create('H:/Threshold/savi_nsdvi.tif', rows, cols, 1,     gdal.GDT_Float32)
outBand = outDs.GetRasterBand(1)

So i want to write all entries in the variable new2 which equal 1 as shapefile.

The white pixels are showing the intersected areas which i want to write to an shapefile.

Afterwards i would like to filter out just those polygons>250m².

  • 1
    I recommend looking at the gdal_sieve.py and gdal_polygonize.py command line tools, they do exactly what you are after. – lpd2 Jan 24 at 10:32
  • Yes thank you i gdal_polygonize.py. I recieve know a polygonized raster. But still with the dimensions of my input raster. How can i manipulate the (rows, cols) that i just recieve values of 1 without zeros. And not in the dimension of the input raster. – RemoteSensing Jan 24 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.