Raster Data Processing¶

Code snippets from CU Earth Data Science Lab Course - Lidar Raster Section

In [1]:

# pyscience imports
import os, sys, glob
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.colors import ListedColormap, BoundaryNorm

import seaborn as sns
plt.style.use("seaborn-darkgrid")
sns.set(style="ticks", context="talk")

import statsmodels.api as sm
import statsmodels.formula.api as smf
# %matplotlib inline
# run for jupyter notebook
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = 'all'

import rasterio as rio
from rasterio.plot import show
from rasterio.plot import plotting_extent
from rasterio.warp import calculate_default_transform, reproject, Resampling
from shapely.geometry import Polygon, mapping
from rasterio.mask import mask

import geopandas as gpd
plt.ion()


# Set standard plot parameters for uniform plotting
plt.rcParams['figure.figsize'] = (8, 8)

# Prettier plotting with seaborn
import seaborn as sns
sns.set(font_scale=1.5)
sns.set_style("white")

Basic Ingestion¶

In [2]:

# Open raster data
lidar_dem = rio.open(
    '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif')

lidar_dem.bounds
#%%
fig, ax = plt.subplots()
show(
    lidar_dem,
    title="Lidar Digital Elevation Model (DEM) \n Boulder Flood 2013",
    ax=ax)
ax.set_axis_off()

Out[2]:

BoundingBox(left=472000.0, bottom=4434000.0, right=476000.0, top=4436000.0)

Out[2]:

<matplotlib.axes._subplots.AxesSubplot at 0x7f24a45350f0>

With context manager¶

In [3]:

#%%
with rio.open(
        '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif'
) as src:
    # convert / read the data into a numpy array: masked= True turns `nodata` values to nan
    lidar_dem_im = src.read(1, masked=True)
    # create a spatial extent object using rio.plot.plotting
    spatial_extent = plotting_extent(src)

print("object shape:", lidar_dem_im.shape)
print("object type:", type(lidar_dem_im))

object shape: (2000, 4000)
object type: <class 'numpy.ma.core.MaskedArray'>

In [4]:

spatial_extent

Out[4]:

(472000.0, 476000.0, 4434000.0, 4436000.0)

Projections / Metadata¶

In [24]:

with rio.open(
        '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif'
) as lidar_dem:
    lidar_dem_im = lidar_dem.read(1, masked=True)
    lidar_dem_im[lidar_dem_im < 0] = np.nan
    lidar_dem.bounds
    myCRS = lidar_dem.crs
    lidar_dem.meta
myCRS

Out[24]:

BoundingBox(left=472000.0, bottom=4434000.0, right=476000.0, top=4436000.0)

Out[24]:

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -3.4028234663852886e+38,
 'width': 4000,
 'height': 2000,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:32613'),
 'transform': Affine(1.0, 0.0, 472000.0,
        0.0, -1.0, 4436000.0)}

Out[24]:

CRS.from_dict(init='epsg:32613')

In [ ]:

Histogram¶

In [6]:

lidar_dem_im.shape
# Remove the `nan` values for plotting
lidar_dem_hist = lidar_dem_im.ravel()
lidar_dem_hist = lidar_dem_hist[~np.isnan(lidar_dem_hist)]

fig, ax = plt.subplots(figsize=(10, 10))
ax.hist(lidar_dem_hist, color='purple')
ax.set(
    xlabel='Elevation (meters)',
    ylabel='Frequency',
    title="Distribution of Lidar DEM Elevation Values")

Out[6]:

(2000, 4000)

Out[6]:

(array([1460307.,  755245.,  615365.,  837086., 1081310., 1135668.,
         687397.,  347037.,  179971.,   58342.]),
 array([1676.21  , 1717.3319, 1758.454 , 1799.5759, 1840.698 , 1881.82  ,
        1922.9419, 1964.064 , 2005.1859, 2046.308 , 2087.43  ],
       dtype=float32),
 <a list of 10 Patch objects>)

Out[6]:

[Text(0, 0.5, 'Frequency'),
 Text(0.5, 0, 'Elevation (meters)'),
 Text(0.5, 1.0, 'Distribution of Lidar DEM Elevation Values')]

Operations on Rasters¶

$CHM = DSM - DTM$

In [7]:

with rio.open(
        '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DSM.tif'
) as lidar_dsm:
    lidar_dsm_im = lidar_dsm.read(1, masked=True)
#%%

# calculate canopy height model
lidar_chm_im = lidar_dsm_im - lidar_dem_im
#%%
# plot the data
fig, ax = plt.subplots(figsize=(10, 6))
chm_plot = ax.imshow(lidar_chm_im, cmap='viridis')
fig.colorbar(chm_plot, fraction=.023, ax=ax)
ax.set_title("Lidar Canopy Height Model (CHM)")
ax.set_axis_off()
#%%

Out[7]:

<matplotlib.colorbar.Colorbar at 0x7f2474362e80>

Out[7]:

Text(0.5, 1.0, 'Lidar Canopy Height Model (CHM)')

In [8]:

dtm_path = '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif'
dsm_path = '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DSM.tif'
with rio.open(dtm_path) as src:
    lidar_dtm_im = src.read(1, masked=True)
    spatial_extent = plotting_extent(src)

with rio.open(dsm_path) as src:
    lidar_dsm_im = src.read(1, masked=True)
    spatial_extent = plotting_extent(src)

lidar_chm_im = lidar_dsm_im - lidar_dtm_im

# View min and max values in the data
print('CHM min value:', lidar_chm_im.min())
print('CHM max value:', lidar_chm_im.max())

CHM min value: 0.0
CHM max value: 26.930054

Exports¶

In [9]:

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

# create dictionary copy
chm_meta = lidar_dem.meta.copy()
# update the nodata value to be an easier to use number
chm_meta.update({'nodata': nodatavalue})
chm_meta

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

out_path = "../spatial/outputs/lidar_chm.tiff"
with rio.open(out_path, 'w', **chm_meta) as outf:
    outf.write(lidar_chm_im_fi, 1)

Out[9]:

(-999.0, -999.0)

Out[9]:

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -999.0,
 'width': 4000,
 'height': 2000,
 'count': 1,
 'crs': CRS.from_dict(init='epsg:32613'),
 'transform': Affine(1.0, 0.0, 472000.0,
        0.0, -1.0, 4436000.0)}

Out[9]:

(2000, 4000)

Binning / Reclassification¶

In [10]:

fig, ax = plt.subplots(figsize=(8, 8))

ax.hist(lidar_chm_im.ravel(), color='purple', edgecolor='white')
ax.set_title(
    "Distribution of Raster Cell Values in the CHM Data", fontsize=16)
ax.set(xlabel="Height (m)", ylabel="Number of Pixels")

Out[10]:

(array([5.846294e+06, 5.079730e+05, 3.710080e+05, 2.406930e+05,
        1.219950e+05, 5.095300e+04, 1.507500e+04, 3.247000e+03,
        4.440000e+02, 4.600000e+01]),
 array([ 0.       ,  2.6930053,  5.3860106,  8.079016 , 10.772021 ,
        13.465027 , 16.158031 , 18.851038 , 21.544043 , 24.23705  ,
        26.930054 ], dtype=float32),
 <a list of 10 Patch objects>)

Out[10]:

Text(0.5, 1.0, 'Distribution of Raster Cell Values in the CHM Data')

Out[10]:

[Text(0, 0.5, 'Number of Pixels'), Text(0.5, 0, 'Height (m)')]

In [11]:

#%% # Histogram
fig, ax = plt.subplots(figsize=(10, 10))
xlim = [0, 25]
ax.hist(
    lidar_chm_im.ravel(),
    color='purple',
    edgecolor='white',
    range=xlim,
    bins=range(*xlim))
ax.set(
    ylabel="Number of Pixels",
    xlabel="Height (m)",
    xlim=[2, 25],
    ylim=[0, 250000])
ax.set_title(
    "Distribution of raster cell values in the DTM difference data\nZoomed in to {} on the x axis"
    .format(xlim),
    fontsize=16)
#%%

Out[11]:

(array([5.448102e+06, 2.445880e+05, 2.201380e+05, 2.061070e+05,
        1.731810e+05, 1.530570e+05, 1.421900e+05, 1.286820e+05,
        1.083210e+05, 8.686700e+04, 6.831900e+04, 5.298100e+04,
        4.023300e+04, 2.987100e+04, 2.082600e+04, 1.369900e+04,
        8.784000e+03, 5.260000e+03, 3.118000e+03, 1.733000e+03,
        8.830000e+02, 4.310000e+02, 2.120000e+02, 9.100000e+01]),
 array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
        17, 18, 19, 20, 21, 22, 23, 24]),
 <a list of 24 Patch objects>)

Out[11]:

[(0, 250000),
 Text(0, 0.5, 'Number of Pixels'),
 (2, 25),
 Text(0.5, 0, 'Height (m)')]

Out[11]:

Text(0.5, 1.0, 'Distribution of raster cell values in the DTM difference data\nZoomed in to [0, 25] on the x axis')

In [12]:

# Define bins that you want, and then classify the data
class_bins = [lidar_chm_im.min(), 2, 7, 12, np.inf]

# You'll classify the original image array, then unravel it again for plotting
lidar_chm_im_class = np.digitize(lidar_chm_im, class_bins)

# Note that you have an extra class in the data (0)
print(np.unique(lidar_chm_im_class))

[0 1 2 3 4]

In [13]:

#%%
# Reassign all values that are classified as 0 to masked (no data value)
# This will prevent pixels that == 0 from being rendered on a map in matplotlib
lidar_chm_class_ma = np.ma.masked_where(
    lidar_chm_im_class == 0, lidar_chm_im_class, copy=True)
lidar_chm_class_ma

# A cleaner seaborn style for raster plots
sns.set_style("white")
# Plot newly classified and masked raster
fig, ax = plt.subplots(figsize=(10, 6))
ax.imshow(lidar_chm_class_ma)
plt.show()
#%%

Out[13]:

masked_array(
  data=[[--, --, --, ..., 1, 1, 1],
        [--, --, --, ..., 1, 1, 1],
        [--, --, --, ..., 1, 1, 1],
        ...,
        [--, --, --, ..., 1, 1, 1],
        [--, --, --, ..., 1, 1, 1],
        [--, --, --, ..., 1, 1, 1]],
  mask=[[ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        ...,
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False],
        [ True,  True,  True, ..., False, False, False]],
  fill_value=999999)

Out[13]:

<matplotlib.image.AxesImage at 0x7f24741b5080>

In [14]:

# Plot data using nicer colors
# Create a list of labels to use for your legend
height_class_labels = [
    "Short trees", "Less short trees", "Medium trees", "Tall trees"
]

colors = ['linen', 'lightgreen', 'darkgreen', 'maroon']

# A path is an object drawn by matplotlib. In this case a patch is a box draw on your legend
# Below you create a unique path or box with a unique color - one for each of the labels above
legend_patches = [
    Patch(color=icolor, label=label)
    for icolor, label in zip(colors, height_class_labels)
]

cmap = ListedColormap(colors)

fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(lidar_chm_class_ma, cmap=cmap)
ax.legend(
    handles=legend_patches,
    facecolor="white",
    edgecolor="white",
    bbox_to_anchor=(1.35, 1))  # Place legend to the RIGHT of the map
ax.set_axis_off()

Out[14]:

<matplotlib.image.AxesImage at 0x7f24741924a8>

Out[14]:

<matplotlib.legend.Legend at 0x7f24741503c8>

Clip / Crop Raster¶

In [15]:

# Open crop extent (your study area extent boundary)
crop_extent = gpd.read_file(
    '../spatial/boulder-leehill-rd/clip-extent.shp')

# Open raster data
chm_path = '../spatial/outputs/lidar_chm.tiff'
lidar_chm = rio.open(chm_path)

#%%
print('crop extent crs: ', crop_extent.crs)
print('lidar crs: ', lidar_chm.crs)
#%% 

crop extent crs:  {'init': 'epsg:32613'}
lidar crs:  EPSG:32613

In [16]:

fig, ax = plt.subplots(figsize=(10, 8))
ax.imshow(
    lidar_chm_im, cmap='terrain', extent=plotting_extent(lidar_chm))
crop_extent.plot(ax=ax, alpha=.8)
ax.set_title("Raster Layer with Shapefile Overlayed")
ax.set_axis_off()

Out[16]:

<matplotlib.image.AxesImage at 0x7f2474118668>

Out[16]:

<matplotlib.axes._subplots.AxesSubplot at 0x7f24741114a8>

Out[16]:

Text(0.5, 1, 'Raster Layer with Shapefile Overlayed')

In [17]:

#%%
# create geojson object from the shapefile imported above
extent_geojson = mapping(crop_extent['geometry'][0])
extent_geojson
#%%
with rio.open(chm_path) as lidar_chm:
    lidar_chm_crop, lidar_chm_crop_affine = mask(
        lidar_chm, [extent_geojson], crop=True)
#%%

Out[17]:

{'type': 'Polygon',
 'coordinates': (((472510.46511627914, 4436000.0),
   (476009.76417479065, 4436000.0),
   (476010.46511627914, 4434000.0),
   (472510.46511627914, 4434000.0),
   (472510.46511627914, 4436000.0)),)}

In [18]:

# Create spatial plotting extent for the cropped layer
lidar_chm_extent = plotting_extent(lidar_chm_crop[0],
                                   lidar_chm_crop_affine)
# Plot your data
fig, ax = plt.subplots(figsize=(10, 8))
ax.imshow(lidar_chm_crop[0], extent=lidar_chm_extent)
ax.set_title("Cropped Raster Dataset")
ax.set_axis_off()
#%%

Out[18]:

<matplotlib.image.AxesImage at 0x7f247408a0b8>

Out[18]:

Text(0.5, 1.0, 'Cropped Raster Dataset')

Reproject¶

In [19]:

lidar_dem = rio.open('../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif')
print(lidar_dem.meta)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -3.4028234663852886e+38, 'width': 4000, 'height': 2000, 'count': 1, 'crs': CRS.from_dict(init='epsg:32613'), 'transform': Affine(1.0, 0.0, 472000.0,
       0.0, -1.0, 4436000.0)}

In [20]:

dst_crs = 'EPSG:3857' # CRS for web meractor 

with rio.open('../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rio.open('../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM_repoject.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rio.band(src, i),
                destination=rio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

In [21]:

def reproject_et(inpath, outpath, new_crs):
    dst_crs = new_crs # CRS for web meractor 

    with rio.open(inpath) as src:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })

        with rio.open(outpath, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rio.band(src, i),
                    destination=rio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest)

In [22]:

crop_extent.crs

Out[22]:

{'init': 'epsg:32613'}

In [23]:

reproject_et( inpath  = '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif', 
              outpath = '../spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM_repoject2.tif', 
              new_crs = crop_extent.crs)