{{ message }}

Instantly share code, notes, and snippets.

# 1kastner/self_drawn_tif_remains_black.py

Last active Apr 6, 2017
 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()

### 1kastner commented Apr 6, 2017

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