gistfile1.txt

#!/usr/bin/env python

from osgeo import gdal, ogr, osr
import numpy as np


def write(  filename,
            rasterOrigin,
            pixelWidth,
            pixelHeight,
            x, y, z,
            array):

    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    options=['COMPRESS=LZW',"BLOCKXSIZE=256", 'TILED=YES']
    memdriver = gdal.GetDriverByName('MEM')
    memRaster = memdriver.Create(filename, cols, rows, 2, gdal.GDT_UInt16 )
    memRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    arrayband = memRaster.GetRasterBand(1)
    arrayband.WriteArray(array)

    zband = memRaster.GetRasterBand(2)
    zs = np.empty_like(array)
    zs.fill(z)
    zband.WriteArray(zs)

    memRaster.BuildOverviews( 'nearest', overviewlist = [2, 4])

    zband.FlushCache()
    arrayband.FlushCache()
    memRaster.FlushCache()

    options=['COMPRESS=LZW',"BLOCKXSIZE=256", 'TILED=YES','COPY_SRC_OVERVIEWS=YES']
    tiffdriver = gdal.GetDriverByName('GTiff')
    tiffds = tiffdriver.CreateCopy(filename, memRaster, options=options )




if __name__ == "__main__":

    # set this to world coordinates
    rasterOrigin = (-123.25745,45.43013)

    # width and height in world coordinates
    pixelWidth = 10
    pixelHeight = 10
    filename = 'test.tif'
    array = np.load('test.npy')

    x = 2
    y = 2

    # a new band with values for Z will be written
    # to the file
    z = 11

    # array needs to be flipped to match TIFF orientation
    reversed_arr = array[::-1]

    write(  filename,
            rasterOrigin,
            pixelWidth,
            pixelHeight,
            x, y, z,
            reversed_arr)