Creating new raster from a polygon and numpy array?

Question

In Python, I have an array of NDVI values for an area and a shapefile that bounds that area. How would you make a new raster file from that?

The NDVI values were created by using the shapefile to clip from a much larger raster and the processed into NDVI, so some info from the original raster is available. However, I would like a raster that only covers the area defined by the boundary, not the original raster.

I am working in Python with GDAL/OGR libraries.

Answer 1

1) Find the extents of the polygon using the ogr library in Python

source_ds = ogr.Open(vector_file)

source_layer = source_ds.GetLayer()

x_min, x_max, y_min, y_max = source_layer.GetExtent()

2) Find the pixle size of the input raster using gdal GetGeoTransform()

dataset = gdal.Open(inputRasterFile)

geotransform = dataset.GetGeoTransform() 

pixelWidth, pixelHeight = geotransform[1], geotranform[5]

3) Write the array into raster using the gdal library

driver = gdal.GetDriverByName('GTiff')

outRaster = driver.Create(newRasterfile, cols, rows, 1, gdal.GDT_Float32)

outRaster.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))

outband = outRaster.GetRasterBand(1)

outband.WriteArray(NDVI_Array)

outRasterSRS = osr.SpatialReference()

outRasterSRS.ImportFromEPSG(4326)

outRaster.SetProjection(outRasterSRS.ExportToWkt())

outband.FlushCache()

If the pixel width and heigth are negative use xmax and xmin instead of xmin and ymin.