subtractraster

from osgeo import gdal, gdalconst
import numpy as np

# open the geoTIFFs from disk
ds1 = gdal.Open('data/file1.tif', gdalconst.GA_ReadOnly)
ds2 = gdal.Open('data/file2.tif', gdalconst.GA_ReadOnly)

# make numpy arrays out of them
arr1 = ds1.ReadAsArray()
arr2 = ds2.ReadAsArray()

# subtract every value from both rasters
result = np.subtract(arr1, arr2)

# get attributes from input raster 1 and set them to the output raster
driver = gdal.GetDriverByName("GTiff")
output = driver.Create('result.tif', ds1.RasterXSize, ds1.RasterYSize, 1, ds1.GetRasterBand(1).DataType)
output.SetGeoTransform(ds1.GetGeoTransform())
output.SetProjection(ds1.GetProjection())

# you can set your desired noData value here if you like
output.GetRasterBand(1).SetNoDataValue(-9999)

# write raster to disk
output.GetRasterBand(1).WriteArray(result)