1

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

closed as unclear what you're asking by BERA, whyzar, lynxlynxlynx, Dan C, MaryBeth Dec 5 '18 at 15:50

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • 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 '18 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 '18 at 15:02
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

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