Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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

This comment has been minimized.

Copy link
Owner Author

@1kastner 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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.