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) outBand.SetNoDataValue(0) outBand.WriteArray(new2)
So i want to write all entries in the variable new2 which equal 1 as shapefile.
Afterwards i would like to filter out just those polygons>250m².