Terrace Long Profiles from Lidar DTM w/ QGIS & Illustrator


Goal of Project

Delineate the flat tops (treads) of glacial terraces along a portion of the Okanogan Valley in north-central Washington. GIS work involves a series of steps that begin with a lidar bare-earth raster (DTM) and ends with a chart showing the longitudinal profiles of the surfaces (Elevation vs. Downstream Distance). This example is for a portion of the Okanogan Valley, WA. The river runs north-south, which makes things much easier than one that runs diagonally. We can use latitude for downstream distance in either case, but its easier on orthogonal valleys.


Software needed

- QGIS 3.26 Buenes Aires (freeware)

- LibreCalc (freeware)

- Adobe Illustrator CC


General set up

- Use the default coordsys set when adding the lidar DTM as first layer (or change as you see fit)

- Output files for this example are saved to "tempgis" folder on my desktop


Steps summary

  1. dtm --> slp (model)

  2. slp -> re (export)

  3. re --> pts (pixels to...)

  4. pts --> elv (sample raster...)

  5. elv --> xy (add x/y...)

  6. manual edits --> RR/RL

  7. xy --> csv (export...)

  8. csv --> chart

  9. chart --> pdf

  10. pdf --> ai


Detailed Steps


1.) Download the DTM tile and HS tile for a portion of the Okanogan Valley from the Washington Lidar Portal

- "okanogan_2015_dtm_12" and "okanogan_2015_dtm_12_hs"

- The DTM raster is for processing. The Hillshade raster is for display only.


The okanogan_2015_dtm_12 tile covers a portion of the Okanogan Valley (Monse-Mallott area).



2.) Open New Project in QGIS

- Store downloaded layers in "TerraceGIS" folder or a folder of your choosing

- Load "okanogan_2015_dtm_12" (DTM) and corresponding "okanogan_2015_dtm_12_hs" (hillshade) into map

- Display both the hillshade

- The default coordsys of data from the Washington Lidar Portal is HARN...something, something...leave it alone unless you have reason to change it


3.) Create a shapefile with several rectangles

- Display Rivermiles points (ECY)

- New Shapefile Layer, "boxes12" in "tempgis" folder, Type = polygon, default coordsys

- Draw a set of rectangles across the DTM, locating tops and bottoms at rivermile intervals (rectangle by extent)

- Right-click to close and enter consecutive number for each beginning with 12 (Ex: 121, 122, 123....1210).

- Save layer when done

Polygons (input boxes) will be used as clippers for the DTM (input raster). Their relatively small size helps reduce processing time. The ideal size for processing is boxes spanning 1-2 km reaches (<2 km downstream distance tiles). I'm using a 2015 MacBook Air laptop (Monterey), which is not fast.



4a.) Option A: Clip DTM to all polygons at once

- This option uses all the box polygons to clip the underlying raster, outputting each to a new temporary layer

- Raster menu > Extraction > Clip Raster by Mask Layer

- Click green recycle button to iterate over this layer

- Match extent of the clipped raster...

- Check box for Keep Resolution of Input Raster

- Click Validate to check

- Do not provide an output file name; the tool will generate the clipped rasters automatically

?? - When processing is complete, export the temp rasters as permanent layers


4b.) Option B: Build a model

- This option creates a model that quickly clips the raster to the box polys and calculates slope for each box...


- Processing Toolbox > Gear button > Create new model > model window opens

- Inputs > Scroll down and drag Vector Features into window

- Description = Box polys

- Geometry type = Polygon...

- OK


Inputs tab...

- Drag Raster Layer into window

- Description = Raster Input Layer


Algorithms tab...

- Drag Clip Raster By Mask Layer from list above into window

- Input = Input Raster Layer

- Mask Layer = Box polys

- No data value = -9999 (click the little up arrow to autofill)

- OK


Algorithms tab...

- Drag Raster Terrain Analysis > Slope into window

- Input = Input Raster Layer

- Elevation Layer = drop down menu > "Clipped (mask)...."

Slope = Slope Layer

- OK


Model Properties...

- Name = Slope for Each Box

- Group = Terrace Tools


Save model > "Slope for Each Box" in "models" folder


...........

Not sure how to add Reclassify tool (reclass slope)...Maybe this?


Algorithms tab...

- Search reclass... > Drag r.reclass into window

- Input raster layer > Algorithm Output > "Slope from algorithm "Slope"...."

- Reclass rules text:

0 thru 3 = 1

3 thru 1000000 = -9999

...........


Run the tool...

- Toolbox > search on Terraces... > "Slope for Each Box"

- Box polys = "Boxes"

- Click the green recycle button to reiterate

- Input raster layer = "okanogan_2015_dtm_12"

- Do not provide a file name (leave as temp)

- Run

- Processing takes a minute

5.) Reclassify slope raster via export (0-3 degrees = flat)

- Select the first newly created slope raster

- Change the data by reclassifying, keeping only pixels with "flat" values, here defined as 0 to 3 degree slope.

- Right click layer > Export > Save As, "121re" in "tempgis" folder

- Set NoData values in window > check box for No Data Values > Hit plus button to add a row > From = 3.01, To = 1000000

- Output file = "121re"


6.) Convert reclassed Slope pixels to points

- Toolbox > Raster Pixels to Points

- Raster layer = "121re"

- Slope values > Change name of filed "VALUE" to "SLOPE_D" (for slope in degrees)

- Output file = "121pts"


Points should look like this if you zoom in. No need to display the layer in the map. Point symbols are shrunk down to 0.2 pt to improve display. Some artifacts (roads) can be seen. Ideally, you would delete road points, but that would take a lot of time. Remember, most roads are located low in the valley, near the modern alluvial floodplain. Glacial terrace treads are much higher, well separated from most roads, reservoirs, and the modern channel of the Okanogan River.



7.) Extract elevation raster values to points

- The tool will append extracted elevation values to the ELEV field in the attribute table for the reclassed slope points shapefile; If this takes a long time, you know there are a lot of flat pixels (points)

- Toolbox > Sample Raster Values

- Input = "121pts"

- Raster layer = "okanogan_2015_dtm_12"

- Output column prefix = ELEV

- Output file = "121elv"


8.) Add X and Y fields to attribute table

- Toolbox > Add X/Y Fields to Layer

- Input = "121elv"

- Processing takes a couple minutes

- Output file = "121xy"


** Steps 5-8 could also be modeled (single layer input, individual outputs) **


Optional: Manually pare down points that are spurious or will not be useful

- There may be obvious groups of points that are not terrace flats (river surface, freeways, lakes, etc.); Or you may wish to analyze points for just one side of the valley (keeping only RR or RL points)

- Make "121xy" points layer active

- Toggle editing to On

- Use select features button

- Delete selected points

- Save Layer Edits button - Rename if you kept river-right or river-left points (Ex: "121xy_RR")

- Toggle editing to Off

- Save project > File menu > Save As....Replace


Optional: Filter elevations below 880'

- Right-click > Filter > ...

- 880' is a value specific to my example data; below 880' is the modern floodplain; I'm reducing processing time by eliminating the modern floodplain; I can represent the river surface in the final compilation later with rivermile elevations


Optional: Remove intermediate files from map

- Layers "121slope", "121re", "121pts", and "121elv" can be removed from map and/or deleted


9.) Export points to CSV file

- CSV, or comma-delimited file, is a generic spreadsheet format useful for moving data between different software; this is the final output file involving QGIS; .csv does not support charts or most formatting; save as another format to keep such things

- Right-click on "121xy" in Table of Contents > Export > Save Features As > Export points to CSV

- File name = "121profile.csv"

- Change export name "SAMPLE_1" to "ELEV"

- Change export name "y" to "LAT"

- Place check marks next to "ELEV" and "LAT"; Don't need "SLOPE_D" (VALUE) or "x" columns

- Uncheck box for Add Saved File to Map (its at the bottom of window)

- Output file = "121profile.csv"

You only need 2 columns to create the profile chart: ELEV and LAT. There should be just 2 check marks on the left.



10.) Save QGIS Project

- Project menu > Save As > change file type to .qgs > Replace...only save a .qgs file (not a .qgz)

- Minimize or close QGIS


11.) Spreadsheet edits

- Delete "121profile.qmd" file

- Open "121profile.csv" in LibreCalc


Format cells...

- Format cells: Elev --> 1 sig fig; Lat --> 5 sig figs


Arrange columns...

- Cut each column > Paste Special > Values and Formats in a new column somewhere to the right

- Arrange the columns so Column A = Lat (x-axis), Column B = Elev (y-axis)


Delete rows (specific to my project)...

- Highlight both columns > Sort on Elev

- Select and delete all rows with elevations 879' and below


Create scatter chart...

- Highlight columns A and B > Chart > Scatter type > Finish


Format data series...

- When chart appears, double-click the points

- Points will turn blue and a formatting window will open

- Symbol = triangle, Color = black, Size = 0.10 proportional, OK

- Click off to unselect Data Series

- Click on the chart to highlight the whole chart > File menu > Export AS PDF


Elevation vs. Downstream Distance plot example.



12.) Adobe Illustrator edits

- Drag PDF into a new AI document

- OK the two pop-up warnings


Remove unneeded elements...

- Use white and black pointers to select and delete invisible bounding box, remove clipping mask, ungroup elements to arrive at isolated points and the various chart elements

- Be careful not to move anything by mistake

- Select chart lines > change to 1pt black line > Group > Shift-slide them to the right out of the way


Slide the chart's lines to the right as shown so you can get at the points. Otherwise, don't mess around with much here. The goal is just to recolor and rescale the triangle symbols (Transform Each) then union the overlapping triangles to form merged clusters (Pathfinder palette), which reduces the number of objects and file size. There are thousands of overlapping triangles you cannot see (useless memory suck).



Format triangles...

- Select all points, but no other elements > No Fill, Grayscale, 100% black, Line weight = 0.25

- File > Save As .ai



Transform each...

- Select only the points > Right-click > Transform > Transform Each > Scale 10% H, 10% V (strokes and effects) > Enter

- If too slow, select some of the points (do it in parts)

- The triangles will now be much smaller, in the same position, and with a line weight of 0.025



Union clusters...

- Select resized clusters > Pathfinder > Add tool to union the triangle clusters; This greatly reduces the number of objects

- File > Save As .ai


14.) Save file

- Group everything

- Save .ai


15.) Compile tiles in Illustrator

- Create a new .ai document for the long profile (compilation)

- Create a regularly-spaced set of ticks for river miles (x-axis)

- Create a vertical scale composed of lines that extend across all river miles; may have to break up into parts

- Open each tile of data; Add each tile to its own layer

- Scale each tile horizontally to match river miles ticks

- Scale each to vertical scale

- Delete extraneous stuff once each tile is set (white arrow)

- Add and format other poster elements

- Prepare for printing on large format plotter (.pdf or .jpg).


Terrace flats appear as discontinuous, gently dipping surfaces. Steeply dipping forms are gullies (ignore). Other random flat spots in the landscape appear as a diffuse scatter of points (ignore). These data are for a 1-mile reach of the Okanogan above 880' elevation. Data from both sides of the valley are shown together. An improvement would be to process separate charts for river-right and river-left and include the modern alluvial plain (points below 880'). For compilation, this single tile would be squished in the horizontal direction, but not the vertical. Red lines are my quick interp of 3 terrace levels.




............................................................................end..........................................................................


Draft notes....


Optional: Separate values into "distance away" columns

- If your river channel runs north-south, then longitude (or easting in UTM) could be used as a measure of "distance away" from the channel. Additional spatial information can be added to the scatter plot. There are other programs (like R) which allow you to do this in spreadsheets easier (i.e., x, y, distance columns). In Numbers, we have to do it with multiple columns that each contain a range of distance values, which corresponds to dot color in the chart.

- Follow the patterns shown here (series and formulas)

- I've used just three distance away columns (series) with a fictitious data set; a few more distance away columns might be good to try; All will be independent series in the chart

- Adjust display colors of points from hot to cold with distance away from river (red = near, blue = far)

Distance away series set up. Values shown are arbitrary. Use decimal degree as the distance unit instead.


Formulas for distance away series columns. Conditional statements are seen at the bottom of each window.

=IF(b2<6,b2,"")

=IF(AND(b2>6,b2<20),b2,"")

=IF(b2>=,b2,"")



Optional: Add a new series for River Mile points

- Download rivermile points

- Select the points for the river

- Add XY



Create label column for the x-axis

- We need to create round values for the labels. Set min x-axis value to 48.24500 and max value to 48.2500. This provides 11 steps between min and max 0.00050 apart. This list will be used to label the axis. Create a column and formula as shown. Your values will differ.

- Create two new columns to the right of the existing data (columns D,E).

- In the next step, we'll point the chart to these values. Follow pattern shown below.


X-axis labels. Create two new columns to the right of the existing data. Type 0.00050 in cell E1. Create formula for cell D2 = D2+$E$2. Drag-copy through cells. Min value for axis is cell D2. Max value is cell D12.



Last 50 Posts
All Posts by Month