self_drawn_tif_remains_black.py

import random

import numpy as np
from matplotlib import pyplot as plt
import rasterio


def draw_map_1():
    # simplified version of what I want to do
    x_size_raster = 500
    y_size_raster = 300
    raster_map = np.ones((x_size_raster, y_size_raster))
    for shape_record in range(200):
        row = random.randint(0, x_size_raster - 1)
        col = random.randint(0, y_size_raster - 1)
        raster_map[row, col] = random.randint(0, 40)
    return raster_map


def draw_map_2():
    # taken from documentation of rasterio
    x_size_raster = 500
    y_size_raster = 300
    x = np.linspace(0.0, 4.0, x_size_raster)
    y = np.linspace(0.0, 3.0, y_size_raster)
    X, Y = np.meshgrid(x, y)
    Z1 = np.exp(-2 * np.log(2) * ((X - 0.5) ** 2 + (Y - 0.5) ** 2) / 1 ** 2)
    Z2 = np.exp(-3 * np.log(2) * ((X + 0.5) ** 2 + (Y + 0.5) ** 2) / 2.5 ** 2)
    raster_map = 10.0 * (Z2 - Z1)
    return raster_map


def draw_map_3():
    # mix map 1 and map 2 and see whether that works
    x_size_raster = 500
    y_size_raster = 300
    x = np.linspace(0.0, 4.0, x_size_raster)
    y = np.linspace(0.0, 3.0, y_size_raster)
    X, Y = np.meshgrid(x, y)
    Z1 = np.exp(-2 * np.log(2) * ((X - 0.5) ** 2 + (Y - 0.5) ** 2) / 1 ** 2)
    Z2 = np.exp(-3 * np.log(2) * ((X + 0.5) ** 2 + (Y + 0.5) ** 2) / 2.5 ** 2)
    raster_map = 10.0 * (Z2 - Z1)
    for shape_record in range(200):
        col = random.randint(0, x_size_raster - 1)
        row = random.randint(0, y_size_raster - 1)
        raster_map[row, col] = random.randint(0, 40)
    return raster_map


def show_map(raster_map):
    print(raster_map.shape)
    print(raster_map.max())
    print(raster_map.min())
    print(raster_map.dtype)
    plt.imshow(np.transpose(raster_map), cmap="hot")
    plt.show()


def save_map(path, raster_map):
    raster_map = np.asarray(raster_map, np.uint8)
    new_dataset = rasterio.open(path, 'w', driver='GTiff',
        height=raster_map.shape[0], width=raster_map.shape[1], count=1,
        dtype=raster_map.dtype, crs='EPSG:4326')
    new_dataset.write(raster_map, 1)
    new_dataset.close()


def open_map(path):
    with rasterio.open(path) as src:
        raster_map = src.read()
        if len(raster_map.shape) == 3:
            raster_map = raster_map[0]  # only use first band
    return raster_map


def demo():
    m1 = draw_map_1()
    m2 = draw_map_2()
    m3 = draw_map_3()
    show_map(m1)
    show_map(m2)
    show_map(m3)
    save_map("m1.tif", m1)
    save_map("m2.tif", m2)
    save_map("m3.tif", m3)
    m1o = open_map("m1.tif")
    m2o = open_map("m2.tif")
    m3o = open_map("m3.tif")
    show_map(m1o)
    show_map(m2o)
    show_map(m3o)

demo()

The issue seems to be with the saved format. QGIS does not like double64-files with a lot of the same color.