Be sure to set your working directory

os.chdir("path-to-you-dir-here/earth-analytics/data")

importrasterioasriofromrasterio.plotimportplotting_extentimportnumpyasnpimportearthpyasetimportmatplotlib.pyplotaspltfrommatplotlib.patchesimportPatchfrommatplotlib.colorsimportListedColormapimportmatplotlib.colorsascolorsimportearthpy.spatialasesimportosimportseabornassnsplt.ion()# Set plot parameters (optional)plt.rcParams['figure.figsize']=(8,8)# prettier plotting with seabornsns.set(font_scale=1.5)sns.set_style("whitegrid")# set working directoryos.chdir(os.path.join(et.io.HOME,'earth-analytics'))

Open and plot the lidar digital elevation model (DEM). Note that when you read the data, you can use the argument masked = True to ensure that the no data values do not plot and are assign nan or nodata.

# open raster datawithrio.open('data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif')aslidar_dem:lidar_dem_im=lidar_dem.read(1,masked=True)# get bounds for plottingbounds=plotting_extent(lidar_dem)

Import Digital Surface Model (DSM)

Next, you will open the digital surface model (DSM). The DSM represents the top of the earth’s surface. Thus, it includes trees, buildings and other objects that sit on the earth.

Canopy Height Model

The canopy height model (CHM) represents the HEIGHT of the trees. This is not an elevation value, rather it’s the height or distance between the ground and the top of the trees (or buildings or whatever object that the lidar system detected and recorded).

Some canopy height models also include buildings, so you need to look closely at your data to make sure it was properly cleaned before assuming it represents all trees!

Calculate difference between two rasters

There are different ways to calculate a CHM. One easy way is to subtract the DEM from the DSM.

DSM - DEM = CHM

This math gives you the residual value or difference between the top of the earth surface and the ground which should be the heights of the trees (and buildings if the data haven’t been “cleaned”).

Data Tip: Note that this method of subtracting 2 rasters to create a CHM may not give you the most accurate results! There are better ways to create CHM’s using the point clouds themselves. However, in this lesson you learn this method as a means to get more familiar with the CHM dataset and to understand how to perform raster calculations in Python.

Before you subtract the two rasters, be sure to check to see if they cover the same area.

# are the bounds the same?print("Is the spatial extent the same?",lidar_dem.bounds==lidar_dsm.bounds),## is the resolution the same ??print("Is the resolution the same?",lidar_dem.res==lidar_dsm.res)

Is the spatial extent the same? True
Is the resolution the same? True

It looks like the bounds and resolution are the same. This means it is safe for you to subtract the two rasters without significant errors or uncertainty introduced.

Export a Raster

You can export a raster file in python using the rasteriowrite() function. Export the canopy height model that you just created to your data folder. You will create a new directory called “outputs” within the colorado-flood directory. This structure allows you to keep things organized, separating your outputs from the data you downloaded.

NOTE: you can use the code below to check for and create an outputs directory. OR, you can create the directory yourself using the finder (MAC) or windows explorer.

Exporting Numpy Arrays to Geotiffs

Next, you need to consider the metdata associated with your chm. Remember that the chm was generated using 2 numpy arrays. Neither of these arrays has spatial data directly associated with it. BUT you do have the rasterio object that has metadata that you can use if you want to assign all of the spatial attributes that are needed to save a usable geotiff file.

You can use the syntax

**dictionary-metadata-object-here

to apply all of the spatial attributes from one of your raster objects, when you write out your new chm raster.

To begin, have a look at the lidar_dem metadata dictionary. Looking at the example below, all of the metadata in that dictionary are the same as what we expect the output chm to have. Thus we can use the metadata as they are.

# fill the masked pixels with a set no data valuenodatavalue=-999.0lidar_chm_im_fi=np.ma.filled(lidar_chm_im,fill_value=nodatavalue)lidar_chm_im_fi.min(),nodatavalue

(-999.0, -999.0)

Then update the metadata dictionary.

# update the metadata to ensure the nodata value is properly documented # create dictionary copychm_meta=lidar_dem.meta.copy()# update the nodata value to be an easier to use numberchm_meta.update({'nodata':nodatavalue})chm_meta

If you want, you can check things like the shape of the numpy array to ensure that it is the same as the width and height of the dem. It should be!

# note the width and height of the dem above. Is the numpy array shape the same?lidar_chm_im_fi.shape

(2000, 4000)

Finally, you can export your raster layer. Below you do the following

you use the same rio.open() syntax that you are used to using except now you specify that you are writing a new file with the ‘w’ argument.

you specify the new file name and destination in the rio.open() function eg: 'data/colorado-flood/spatial/outputs/lidar_chm.tiff'

you specify the metadata as an “unpacked” dictionary using **lidar_dem.meta - doing this allows you to NOT have to specify EACH and EVERY metadata element individually in your output statement - which would be tedious!

finally you write the file. output_file.write(your-object-name, layer) Notice that when you make this call you specify both the object name and the layer that you wish to write to a new file. Also notice the outf is the name of the rasterio object as defined below.