from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg


# import data
ds = gdal.Open('data/dem.tif')
gt = ds.GetGeoTransform()
proj = ds.GetProjection()

band = ds.GetRasterBand(1)
array = band.ReadAsArray()
plt.figure()
plt.imshow(array)
plt.show()


# manipulate
binmask = np.where((array >= np.mean(array)),1,0)
plt.figure()
plt.imshow(binmask)
plt.show()

driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("binmask.tiff", xsize = binmask.shape[1], ysize = binmask.shape[0], bands = 1, eType = gdal.GDT_Int16)
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(binmask)
outband.SetNoDataValue(np.nan)
outband.FlushCache()

outband = None
outds = None
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg


# import data
ds = gdal.Open('data/dem.tif')
gt = ds.GetGeoTransform()
proj = ds.GetProjection()

band = ds.GetRasterBand(1)
array = band.ReadAsArray()
plt.figure()
plt.imshow(array)
plt.show()


# manipulate
binmask = np.where((array >= np.mean(array)),1,0)
plt.figure()
plt.imshow(binmask)
plt.show()

driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("binmask.tiff", xsize = binmask.shape[1], ysize = binmask.shape[0], bands = 1, eType = gdal.GDT_Int16)
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(binmask)
outband.SetNoDataValue(np.nan)
outband.FlushCache()

outband = None
outds = None
# import data # manipulate