Question 1¶
import rioxarray as rxr
import xarray as xr
import geopandas as gpd
import earthpy.plot as ep
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.colors import ListedColormap, BoundaryNorm
import earthpy as et
import numpy as np
import seaborn as sns
sns.set(font_scale=1.5, style="whitegrid")
Importing necessary libraries for project.
myTif = rxr.open_rasterio('../earth-analytics/data/colorado-flood/spatial/outputs/pre-flood-chm.tif',
masked=True).squeeze()
myTif
Opening pre-flood canopy height model and viewing data array.
print("The CRS for this data is:", myTif.rio.crs)
print("The no data value is:", myTif.rio.nodata)
print("The spatial extent is:", myTif.rio.bounds())
The CRS is correct.
myTif.plot(cmap=plt.cm.BuPu_r)
plt.show()
Final plot with blue-purple color bar.
Question 2A¶
pre_DTM = rxr.open_rasterio('../earth-analytics/data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DTM.tif',
masked=True).squeeze()
pre_DSM = rxr.open_rasterio('../earth-analytics/data/colorado-flood/spatial/boulder-leehill-rd/pre-flood/lidar/pre_DSM.tif',
masked=True).squeeze()
post_DTM = rxr.open_rasterio('../earth-analytics/data/colorado-flood/spatial/boulder-leehill-rd/post-flood/lidar/post_DTM.tif',
masked=True).squeeze()
post_DSM = rxr.open_rasterio('../earth-analytics/data/colorado-flood/spatial/boulder-leehill-rd/post-flood/lidar/post_DSM.tif',
masked=True).squeeze()
This cell loads the terrain and surface models from before and after the flood.
# canopy height models (CHM)
# pre
pre_CHM = pre_DSM - pre_DTM
#post
post_CHM = post_DSM - post_DTM
This cell creates two canopy height models: one from before the flood and one from after.
canopy_difference = post_CHM - pre_CHM
This cell creates a data array of the canopy height difference by subtracting the height of the canopy before the flood from the height of the canopy after the flood.
print('min value:', np.nanmin(canopy_difference))
print('max value:', np.nanmax(canopy_difference))
Here I am checking the min and max values of the data which I will need to know when creating the bins in the next step.
class_bins = [-25, 0, 25]
canopy_difference_class = xr.apply_ufunc(np.digitize,
canopy_difference,
class_bins)
In this cell I am creating a classified version of the data. The value '0' will be included in the category 'Positive Change'.
canopy_difference_class_ma = canopy_difference_class.where(canopy_difference_class != 3)
In this cell, I am creating a new data array object by selecting all data except for the data classified as'nodata'.
class_labels = ["Negative Change",
"Positive Change"]
colors = ['maroon',
'gold']
cmap = ListedColormap(colors)
class_bins = [.5, 1.5, 2.5]
norm = BoundaryNorm(class_bins,
len(colors))
# Plot newly classified and masked raster
f, ax = plt.subplots(figsize=(10, 5))
im = canopy_difference_class_ma.plot.imshow(cmap=cmap,
norm=norm,
# Turn off colorbar
add_colorbar=False)
# Add legend using earthpy
ep.draw_legend(im,
titles=class_labels)
ax.set(title="Classified Lidar Pre/Post-Flood Canopy Height Difference Model \n Derived from NEON AOP Data")
ax.set_axis_off()
plt.show()
In this cell I set the class labels, colors and bins for plotting the classified and masked data. Then I plotted it. The reason I included 0 height difference in the 'Positive Change' category is because I think that this would be considered a positive outcome in a real world situation. I also added a legend and a title, and turned off the axis.
Question 2B¶
terrain_change = post_DTM - pre_DTM
This cell creates a data array of the terrain height difference by subtracting the height of the terrain before the flood from the height of the terrain after the flood.
print('min value:', np.nanmin(terrain_change))
print('max value:', np.nanmax(terrain_change))
Here I am checking the min and max values of the data which I will need to know when creating the bins in the next step.
class_bins = [-20, 0, 50]
terrain_change_class = xr.apply_ufunc(np.digitize,
terrain_change,
class_bins)
In this cell I am creating a classified version of the data. The value '0' will be included in the category 'Positive Change'.
terrain_change_class_ma = terrain_change_class.where(terrain_change_class != 3)
In this cell, I am creating a new data array object by selecting all data except for the data classified as 'nodata'.
class_labels = ["Negative Change",
"Positive Change"]
colors = ['maroon',
'gold']
cmap = ListedColormap(colors)
class_bins = [.5, 1.5, 2.5]
norm = BoundaryNorm(class_bins,
len(colors))
# Plot newly classified and masked raster
f, ax = plt.subplots(figsize=(10, 5))
im = terrain_change_class_ma.plot.imshow(cmap=cmap,
norm=norm,
# Turn off colorbar
add_colorbar=False)
# Add legend using earthpy
ep.draw_legend(im,
titles=class_labels)
ax.set(title="Classified Lidar Pre/Post-Flood Terrain Change Model \n Derived from NEON AOP Data")
ax.set_axis_off()
plt.show()
In this cell I set the class labels, colors and bins for plotting the classified and masked data. Then I plotted it. Again, the reason I included 0 height difference in the 'Positive Change' category is because I think that this would be considered a positive outcome in a real world situation. I also added a legend and a title, and turned off the axis.