2

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.

  • You don't state which GIS system you are using? If you have access to ArcMap then there is an arcpy function NumPyArrayToRaster which you could use. – Hornbydd Aug 26 '16 at 8:49
2

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.

| improve this answer | |

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.