2

I would like to create a binary mask using the height values from a Digital Elevation Model, and export this out as a tif. For example, generate a mask for areas over 50m

My code so far:

from osgeo import gdal
import numpy as np
ds = gdal.Open("DEM.tif")
DEM = np.array(ds.GetRasterBand(1).ReadAsArray())

print(DEM.shape)
print(DEM.size)
print(DEM)

This results in reading the raster as an Array

The next step is for me to create a binary mask of values over, say, 1000m, and export this as a tif. file

2
  • 1
    What have you tried so far? Have you looked at any tutorials and are stuck at a specific point of the workflow?
    – Erik
    Dec 5, 2018 at 14:07
  • @Erik, please see my edited post. I have read the raster file as a Numpy Array. I would now like to create a binary mask of the areas over 1000m
    – Rhodri_22
    Dec 5, 2018 at 15:02

1 Answer 1

3

You can create a raster mask using this:

from osgeo import gdal
import numpy as np
ds = gdal.Open("Norway_WD4O_Copy.tif")
ulx, xres, xskew, uly, yskew, yres =ds.GetGeoTransform()
band_ds = ds.GetRasterBand(1)
DEM = band_ds.ReadAsArray()

Mask_Zero = DEM < 1000 #In units of the raster, e.g. meters
Ones_array = np.ones((DEM.shape[0],DEM.shape[1]),dtype=np.float32)

Ones_array[Mask_Zero] = 0
#This will create a mask of 1 where the DEM is equal or above 1000m and 0 when below.

out_raster = gdal.GetDriverByName("GTIFF").Create(r"path\to\folder\name.tif",DEM.shape[1],DEM.shape[0],1,gdal.GDT_Float32)
out_raster.SetGeoTransform((ulx, xres, xskew, uly, yskew, yres))
band = out_raster.GetRasterBand(1)
band.WriteArray(Ones_array) 
out_raster = None
band = None

And if you want it as a vector you can use gdal polygonize

0

Not the answer you're looking for? Browse other questions tagged or ask your own question.