Ignitions: Lines following contours

zacharycope0 commented 2 years ago

Feature Description When there is topography the ignitors with consider both wind and aspect when putting down lines. As Alex showed in the meeting he has code to build ignitions along contours, but he hasn't automated a way to build the contours. I would like to make a function that builds contour lines from topo.dat. I would like to lines to be numbered from 1-n starting with the lines at a higher elevation.

Getting Started First work on building and numbering the lines.

Then reach out to Alex to see how it fits into his workflow.

Screenshots Add any other context or screenshots about the feature request here.

mathfire commented 2 years ago

I was messing around with building from a DEM raster today and had some good luck. I am trying not to use a crazy number of packages to get it to work. Here are some highlights:

This routine relies on matplotlib and rasterio.
I am relying on rasterio to import a .tif Digital Elevation Model (not a topo.dat... yet)
ct_gap_m is the contour spacing in meters.
I used matplotlib.plot.contour() to return the contour lines,
Then I unpack the contours with a loop.

Here's the code snippet:

# Episode 1
# Working towards seamless rendering of Digital Elevation Model
# Also obtaining contours without use of Arc or QGIS products

# Copy-and-paste below to run in Python console
# exec(open('C:\\Users\\FireScience\\Documents\\QUIC-fire\\ignPad\\dem_read_ep1.py').read())

import os
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
from rasterio.plot import show_hist

fp = r"C:\\Users\\FireScience\\Documents\\QUIC-Fire\\ignPad\\GISdata\\Konza_SHP\\DEM_Konza.tif"

raster = rasterio.open(fp)

# Vital statistics about raster
print('# of bands = ' + str(raster.count))
print('# of rows = ' + str(raster.height))
print('# of columns = ' + str(raster.width))
print('bbox listed below')
print('-----------------')
print(raster.bounds)
print('LL (Lower Left) corner listed below')
print('-----------------------------------')
print(raster.transform * (raster.width, raster.height))

# Matplotlib plotting
fig, axs = plt.subplots(1, 2, figsize=(16, 9))
show(raster.read(1), cmap='Greys_r', transform=raster.transform, ax=axs[0])
#show(raster, ax=axs[0])
# Histogram
show_hist(raster, bins=50, lw=0.0, stacked=False, alpha=0.3,
          histtype='stepfilled', title="Histogram", ax=axs[1])
plt.show()

# Work on developing contours
ct_gap_m = 10.0
matrix = raster.read(1)
#matrix = ((np.floor(matrix + ct_gap_m + 1.0) / ct_gap_m) * ct_gap_m) - ct_gap_m
num_intervals = int((matrix.max() - matrix.min()) / ct_gap_m)
fig, axs = plt.subplots(1, 2, figsize=(16, 9))
ct = axs[0].contour(matrix, np.linspace(matrix.min(), matrix.max(), num_intervals))
contours = []
for cc in ct.collections:
    paths = []
    for pp in cc.get_paths():
        xy = []
        for vv in pp.iter_segments():
            xy.append(vv[0])
        paths.append(np.vstack(xy))
    contours.append(paths)
for line in contours:
    if line:
        line = line[0]
        axs[1].plot(line[:, 0], line[:, 1])
        print('contour length = {0}'.format(len(line)))
    else:
        print('found an empty contour!')
plt.show()

cbonesteel commented 2 years ago

@mathfire I finally got around to looking at this. This does look promising! I haven't exactly messed around with rasterio yet so some guidance on how to get a .tif file that would be able to run in this code would be appreciated. I'm hoping that tomorrow I can take a look into this code a bit more and see how it can be implemented into the standardized workflow.

mathfire commented 2 years ago

@cbonesteel Great! One useful online clearinghouse for DEM data is: https://apps.nationalmap.gov/downloader/ This is the so-called "The National Map" offered by the USGS. To use it to find those DEMs (.tif and .tiff) involves defining a bounding box/extent, then searching for "Elevation Products (3DEP)" and download. For my case, I am finding little sections of terrain for the ignitions that Daniel asks me to draw and a list of escaped RX fires of interest.

Otherwise, I realized the code snippet I shared had poor commenting. So sorry about that! I like to lead by example so I def fell down there. Realistically, n lines of code should be at least 2n lines of written work. Each line with a comment. Don't know if that is how you practice or were taught, but it makes collaboration so much easier. Here is a better formatted snippet:

# Episode 1
# Working towards seamless rendering of Digital Elevation Model
# Also obtaining contours without use of Arc or QGIS products

# Copy-and-paste below to run in Python console
# exec(open('/home/aleks/PycharmProjects/pythonProject/dem_read_ep1.py').read())

# Supporting libraries and dependencies
import os
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
from rasterio.plot import show_hist

# Create a read file pointer to a DEM (Digital Elevation Model)
fp = r"/home/aleks/PycharmProjects/pythonProject/Konza_SHP/DEM_Konza.tif"
# Stream in the raster grid from the file
raster = rasterio.open(fp)
# Report to the console some vital statistics about raster
print('# of bands = ' + str(raster.count))      # Number of color bands
print('# of rows = ' + str(raster.height))      # Height (# of rows of raster)
print('# of columns = ' + str(raster.width))    # Width (# of columns of raster)
print('bbox listed below')
print('-----------------')
print(raster.bounds)                            # A header, then the bounding box
print('LL (Lower Left) corner listed below')
print('-----------------------------------')
print(raster.transform * (raster.width, raster.height))     # A header, then the LL corner

# Matplotlib plotting
fig, axs = plt.subplots(1, 2, figsize=(16, 9))
# Show the raster on the left subplot
show(raster.read(1), cmap='Greys_r', transform=raster.transform, ax=axs[0])
#show(raster, ax=axs[0])        # Alternate command to so the same
# Histogram of the colors
show_hist(raster, bins=50, lw=0.0, stacked=False, alpha=0.3,
          histtype='stepfilled', title="Histogram", ax=axs[1])
# Show current figure
plt.show()

# Establish the contour spacing in meters
ct_gap_m = 10.0
# Save a copy of the raster
matrix = raster.read(1)
# "Integerize" the DEM for evenly spaced level sets
#matrix = ((np.floor(matrix + ct_gap_m + 1.0) / ct_gap_m) * ct_gap_m) - ct_gap_m
# Calculate the number of contour intervals
num_intervals = int((matrix.max() - matrix.min()) / ct_gap_m)
# Next figure to draw
fig, axs = plt.subplots(1, 2, figsize=(16, 9))
# Save to a local variable and draw the contours
ct = axs[0].contour(matrix, np.linspace(matrix.min(), matrix.max(), num_intervals))

# Loop to extract individual contour lines
# Empty container where the contour line coordinates will be stored
contours = []
# Loop through each member object of the ct return from Matplotlib contour
for cc in ct.collections:
    # Empty container for paths within the cc member objects
    paths = []
    # Loop through the paths
    for pp in cc.get_paths():
        # Initialize space for coordinates
        xy = []
        # Loop through segments of the paths
        for vv in pp.iter_segments():
            # Add coordinates to the list
            xy.append(vv[0])
        # Stack the coordinates and add them to this path
        paths.append(np.vstack(xy))
    # Add the path (with coordinates) to the contour
    contours.append(paths)
# For each contour line
for line in contours:
    # If there are coordinates
    if line:
        # Parse the coordinates
        line = line[0]
        # PLot the coordinates
        axs[1].plot(line[:, 0], line[:, 1])
        # Report the length of this contour
        print('contour length = {0}'.format(len(line)))
    # In absence of contours...
    else:
        # ... report an empty
        print('found an empty contour!')
# Show current figure
plt.show()

QUIC-Fire-TT / ttrs_quicfire

Ignitions: Lines following contours #6