Using raster calculator with GDAL in python?

Question

i am using an hyperspectral image with 158 bands. I want to calculate raster. i am using the following python script

import gdal # Import GDAL library bindings
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
import pylab as plt
import numpy as np
import xlrd
# The file that we shall be using
# Needs to be on current directory
filename = ('C:/Users/KIFF/Desktop/These/data/Hyperion/10th_bandmathref')
outFile = ('C:/Users/KIFF/Desktop/These/data/Hyperion/Math')
XLS=('C:/Users/KIFF/Desktop/These/data/Coef/bcoef.xlsx')
wb = xlrd.open_workbook(XLS)
sheet = wb.sheet_by_index(0)
sheet.cell_value(0, 0)
g = gdal.Open(filename, GA_ReadOnly)
# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
 print ("Problem opening file %s!" % filename)
else:
 print ("File %s opened fine" % filename )
#band_array = g.ReadAsArray()
#print(band_array)
print ("[ RASTER BAND COUNT ]: ", g.RasterCount)
for band in range( g.RasterCount ):
 print (band)
 band += 1
 outFile = ('C:/Users/KIFF/Desktop/These/data/Results/Temp/Math_1_sur_value'+str(band)+'.tiff')
 #print ("[ GETTING BAND ]: ", band )
 srcband = g.GetRasterBand(band)
 if srcband is None:
 continue
 data1 = BandReadAsArray(srcband).astype(np.float)
 print(data1)
 # for i in range(3,sheet.nrows):
 b=sheet.cell_value(band+2,1)
 #print(b)
 dataOut = (1/data1)
 driver = gdal.GetDriverByName("ENVI")
 dsOut = driver.Create(outFile, g.RasterXSize, g.RasterYSize, 1)
 CopyDatasetInfo(g,dsOut)
 bandOut=dsOut.GetRasterBand(1)
 BandWriteArray(bandOut, dataOut)

for the print(data1) i got only some "1", but the real values are some floats

0
[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
1
[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
2

Pixel value 0,139200

Can you find the mistake?

Answer 1

Try setting the raster data type - such as gdal.GDT_Float32

dsOut = driver.Create(outFile, g.RasterXSize, g.RasterYSize, 1, gdal.GDT_Float32)

• Kais Tounsi: Thanks. It's working!!!!