# Принимаем, что у нас уже есть все матрицы lats, lons, div, ...
# Покажем, как через GCPs сохранить матрицу div в geotiff

import gdal, osr

# Готовим GCPs

gcps = []

for i in range (0,div.shape[0]):
    for j in range(0,div.shape[1]):
        gcps.append(gdal.GCP(lons[i,j], lats[i,j], 0, i, j))

# Определяем систему координат GCPS как WGS84
gcp_srs = osr.SpatialReference()
gcp_srs.ImportFromEPSG(4326)
gcp_wkt = gcp_srs.ExportToWkt()

# Создаём временный gdal-датасет в памяти
driver = gdal.GetDriverByName("MEM")
dataType = gdal.GDT_Float32
dataset = driver.Create('', div.shape[1], div.shape[0], 1, dataType)

# Пишем в него нужную матрицу
dataset.GetRasterBand(1).WriteArray(div)

# Устанавлиаем GCPs
dataset.SetGCPs(gcps, gcp_wkt)

# Проецируем в нормальный вид (переход от GCP к проецированному растру)
gdal.Warp('output.tif',dataset)

# Можно сразу в этой же команде задать проекцию, другие размеры пикселя и т.д., напр.
# gdal.Warp('output.tif',dataset,xRes=...,yRes=...,dstSRS=...)

del dataset