from osgeo import gdal
import numpy

g = gdal.Open ("baikal_subset.tif")
if g is None:
    raise IOError, "Couldn't open baikal_subset.tif"

b3 = g.GetRasterBand(3).ReadAsArray().astype(np.float32)
b4 = g.GetRasterBand(4).ReadAsArray().astype(np.float32)

ndvi = (b4 - b3)/(b4 + b3)

drv = gdal.GetDriverByName ( "GTiff" )
dst_ds = drv.Create ( "output_file.tif", g.RasterXSize, g.RasterYSize, 1, 
       gdal.GDT_Float32, options=["COMPRESS=LZW"] )
dst_ds.GetRasterBand(1).WriteArray ( dst_ds.astype (np.float32) )
dst_ds = None