Analysing along a feature: oblique frames

A fault, a trunk river, a coastline or a mountain range rarely runs north–south or east–west. Analysed on a UTM grid, a long diagonal strip sits in a bounding box that is mostly not the strip: a 100 km × 10 km strip at 45° fills a 78 km × 78 km box, six times the cells. Orientations come out relative to north, when what you usually want is relative to the feature.

lsdtt3 can instead put the DEM in an oblique frame: a grid whose x axis runs along the feature. The DEM is still an ordinary north-up GeoTIFF; only its CRS is different (a Hotine Oblique Mercator projection centred on the feature), so every lsdtt3 tool reads and analyses it unchanged. When the analysis is done, products such as basins, channel masks and slope are rotated back onto a UTM grid. How the frame is built, and how well results agree between frames and UTM, is in How oblique frames work at the end of this page.

This page works through two examples with the command-line tools:

a strip along the San Andreas Fault in the Carrizo Plain, California, fetched straight into its frame, analysed, rotated back and plotted;
the Wind River Range, Wyoming: an existing UTM DEM rotated into the frame of a polygon drawn around the range, then profiled along its axis.

The same steps from Python are in Oblique frames in Python.

What a frame looks like

+x runs along the feature and the frame centre is (0, 0). The axis azimuth is stored folded into (−90°, 90°], so for a feature trending north-west to south-east (azimuth −45°) +x points north-west and the map has its north-west end on the right.
Up is not north. In that −45° frame, grid up points south-west, and true north is 135° clockwise from grid up. Tools that write bearings print this at start-up:

Grid north is 135.4 deg anticlockwise of true north at the grid centre;
bearings are relative to grid north (bearings_relative_to_true_north is false).

Bearings are relative to the grid unless you ask otherwise: aspect, channel-link bearings, basin orientation and the other bearing outputs are clockwise from grid up, so 90° means "along the feature". Set bearings_relative_to_true_north: true in any tool that writes bearings to get them from true north instead; each value is converted at its own location. Angles between two directions (junction angles, deflections) are the same either way, and hillshade illumination always stays relative to the grid (315° lights the map from its top left).
A rotation frame file is written next to every DEM put into a frame (<name>_rotation_frame.txt). It records the frame and the UTM grid to rotate products back onto. Keep it with the DEM: rotate-back reads it.

Example 1: a fault strip

The feature file

The fault is a line in any vector format and CRS. Here it is an approximate trace of the San Andreas Fault across the Carrizo Plain, digitised for this example, saved as carrizo_fault.geojson:

{"type": "FeatureCollection", "name": "carrizo_fault",
 "crs": {"type": "name", "properties": {"name": "EPSG:4326"}},
 "features": [{"type": "Feature", "properties": {"name": "San Andreas Fault, Carrizo Plain (approximate)"},
   "geometry": {"type": "LineString", "coordinates": [[-120.02, 35.40], [-119.80, 35.22], [-119.55, 35.02]]}}]}

A line (or a multi-part line) is framed by its two furthest end points; a file with two points uses them as the axis ends; a polygon is framed by its minimum-area rectangle, with the longest side as the axis (Example 2).

1. Fetch the strip into its frame

# fetch.param: an 8 km wide strip along the fault, straight into its frame.
fetch_from_linear_feature: true
fetch_feature_fname: carrizo_fault.geojson
linear_feature_fetch_width: 8000
grid_spacing: 30
source: cop30_aws
write_fname: carrizo

lsdtt-fetch-raster fetch.param

Centre  : -119.7844525, 35.21023104
Azimuth : -45.42894765 deg
Size    : 60720 m x 8040 m (2024 x 268 cells)
Wrote: ./carrizo_oblique_DEM.tif

The source tiles are warped straight into the frame, so the DEM is resampled once. The strip is the line's length plus an end buffer (linear_feature_fetch_end_buffer, default ten cells at each end) by the width. Next to the DEM is carrizo_oblique_DEM_rotation_frame.txt, whose rotate-back target is UTM zone 11N (EPSG:32611), the zone of the frame centre, on the 30 m lattice a plain fetch uses. rotate_after_fetch: false would instead fetch the strip's bounding box on that UTM grid.

The frame holds 2024 × 268 = 542,432 cells. The same strip rotated back onto UTM 11N needs a 1597 × 1648 grid, 2,631,856 cells.

2. Analyse in the frame

Any tool runs on the frame DEM as on any projected DEM. Here, complete second-order basins (in an 8 km strip, third-order basins were all cut by the strip's edges), and slope and aspect:

# basins.param
read_fname: carrizo_oblique_DEM.tif
write_fname: carrizo
threshold_contributing_pixels: 500
select_all_basins_of_stream_order: true
selected_stream_order: 2
filter_incomplete_basins: true
write_basins_raster: true
write_basin_outlines_as_polygon: true

# metrics.param: aspect is clockwise from grid up, 90 = along the fault.
read_fname: carrizo_oblique_DEM.tif
write_fname: carrizo
polyfit_metrics: slope,aspect_degrees
polyfit_lambda_metres: 90
write_hillshade_raster: true

lsdtt-basin-extraction basins.param      # 65 basins: carrizo_junction_basins.tif + _polygons.fgb
lsdtt-surface-metrics metrics.param      # carrizo_direct_polyfit_slope.tif, _aspect_degrees.tif

3. Rotate the products back

Rotate-back is a separate mode of lsdtt-raster-preprocessing. List the rasters one per line, with near for anything categorical (basin IDs, masks, junction numbers, lithology); the default method is bilinear:

# rotate_back_list.txt
carrizo_junction_basins.tif, near
carrizo_direct_polyfit_slope.tif

# rotate_back.param
run_rotate_back: true
rotation_frame_fname: carrizo_oblique_DEM_rotation_frame.txt
rotate_back_file_list: rotate_back_list.txt
write_fname: carrizo

lsdtt-raster-preprocessing rotate_back.param

Each raster is written as <name>_rotated_back.tif on EPSG:32611 at 30 m, here 1597 × 1648 cells, with all 65 basin IDs. Vector outputs such as carrizo_junction_basin_polygons.fgb need no rotating back: they carry their CRS, and QGIS, geopandas and lsdviztools3 reproject them.

Rotate back products, never the DEM. Rotating a DEM into a frame and back resamples it twice: for the 30 m Copernicus DEM of Example 2 (cubic both ways, 5.6 million cells), the round trip changed elevations by 0.36 m RMS, 1.3 m at the 99th percentile and up to 18 m on cliffs. Keep the original DEM for anything that needs elevations on the UTM grid.

4. Figures

lsdviztools3 draws a frame DEM as a rectangle in frame km, with axes labelled along and across the feature, a title giving the frame azimuth, the UTM grid drawn over it as dashed lines with their eastings and northings, and a north arrow that points to true north at its own position:

from lsdviztools3.io.raster import get_crs, load_raster, to_coordinates, to_gmt_grid
from lsdviztools3.render.base import grid_to_display
from lsdviztools3.render.decorations import add_north_arrow
from lsdviztools3.workflows import hillshade_with_vectors


def map_with_arrow(base, layers, png, arrow_size=0.07, **kwargs):
    """A base raster + vector layers + a north arrow pointing to true north."""
    fig = hillshade_with_vectors(base, layers, legend=False, **kwargs)
    grid = to_gmt_grid(to_coordinates(load_raster(base), "local"))
    crs = get_crs(grid)
    _, extent = grid_to_display(grid, crs)
    add_north_arrow(fig, extent, crs, anchor="BR", label_halo="white",
                    size_frac=arrow_size)
    fig.savefig(png, dpi=200)


fault = {"path": "carrizo_fault.geojson", "kind": "line", "color": "red", "width": "1.5p"}
basins = {"path": "carrizo_junction_basin_polygons.fgb", "kind": "polygon",
          "color": "blue", "width": "0.6p"}

# The strip in its frame: the fault runs along +x.
map_with_arrow("carrizo_oblique_DEM.tif", [basins, fault], "carrizo_frame_map.png",
               arrow_size=0.25)

# Slope rotated back onto UTM 11N; the basin polygons, still in the frame CRS,
# are reprojected onto the map.
map_with_arrow("carrizo_direct_polyfit_slope_rotated_back.tif", [basins, fault],
               "carrizo_utm_map.png", shade=False, base_cmap="viridis")

The graticule is on by default for frame maps and off for all others; MapStyle(graticule=False) turns it off, and graticule_interval_km and reference_crs set its spacing and CRS. The north arrow is turned by the grid convergence on every projected map: by up to about 3° on UTM, by 135° in this frame.

Example 2: an elongated range from an existing DEM

When the DEM is already on disk, lsdtt-raster-preprocessing rotates it into a frame in place of its usual reprojection step. Here the DEM is a plain UTM fetch of the Wind River Range, and the frame comes from a polygon drawn around the range (124 km × 40 km, saved as wind_river_range.geojson; any vector format and CRS):

{"type": "FeatureCollection", "name": "wind_river_range",
 "crs": {"type": "name", "properties": {"name": "EPSG:4326"}},
 "features": [{"type": "Feature", "properties": {"name": "Wind River Range (approximate)"},
   "geometry": {"type": "Polygon", "coordinates": [[[-109.751, 43.507], [-110.148, 43.292],
     [-109.245, 42.392], [-108.854, 42.607], [-109.751, 43.507]]]}}]}

1. The UTM DEM

# fetch.param: an ordinary UTM DEM covering the polygon's bounding box.
bounds_fname: wind_river_range.geojson
grid_spacing: 30
source: cop30_aws
write_fname: windriver

lsdtt-fetch-raster fetch.param    # windriver_DEM.tif: 3593 x 4191 cells, EPSG:32612

2. Rotate it into the range's frame

# rotate.param
read_fname: windriver_DEM.tif
write_fname: windriver
rotate_to_linear_feature: true
rotation_feature_fname: wind_river_range.geojson
write_fill_raster: true
write_hillshade_raster: true

lsdtt-raster-preprocessing rotate.param

Centre  : -109.4957472, 42.94983547
Azimuth : -36.30276092 deg
Size    : 124920 m x 40680 m (4164 x 1356 cells of 30 m)
Rotate back : EPSG:32612 at 30, lattice origin (568830, 4693650)

The rotated DEM is windriver_oblique.tif (4164 × 1356 = 5,646,384 cells, against 15,058,263 for the UTM DEM), and every later step in the same run (here the fill and the hillshade) uses it. Because the source DEM is projected in metres, the frame keeps its 30 m cells and the rotate-back target is its own grid: products rotated back land on the source DEM's pixels. The frame is on the source DEM's own datum, so a NAD83 DEM gets a NAD83 frame.

Other ways to define the frame:

rotate_to_linear_feature with a line needs linear_feature_rotation_width (metres); a polygon takes its width from the polygon.
rotate_to_data_footprint: true frames the DEM's valid cells (the minimum-area rectangle around them) instead of a feature file, for a DEM already clipped to an elongated area. A footprint with no clear long axis is warned.
linear_feature_rotation_end_buffer adds length beyond the axis ends (all four sides for a polygon or a footprint); the default, −1, is ten cells.

3. A swath along the range

In the frame the range axis is the x axis, so a swath baseline along the range is the line y = 0. lsdtt-swath-profile reads the baseline's coordinates as frame metres:

{"type": "FeatureCollection", "name": "range_axis",
 "features": [{"type": "Feature", "properties": {},
   "geometry": {"type": "LineString", "coordinates": [[-62000.0, 0.0], [62000.0, 0.0]]}}]}

# swath.param
read_fname: windriver_oblique.tif
write_fname: windriver
baseline_vector_fname: range_axis.geojson
half_width_metres: 15000
bin_width_metres: 2000
write_distance_rasters: false

lsdtt-swath-profile swath.param    # windriver_swath_profile.csv

and plot it, with the map of the range in its frame:

from lsdviztools3.charts import plot_swath
from lsdviztools3.io.swath import load_swath

fig, ax = plot_swath(load_swath("windriver_swath_profile.csv"))
ax.set_xlabel("Distance along the range axis (km)")
fig.savefig("windriver_swath.png", dpi=200, bbox_inches="tight")

Rules of thumb

Analyse in the frame; bring back products, never the DEM.
Bearings are from grid up unless bearings_relative_to_true_north: true.
Do not set force_reproject on a frame DEM: it would warp it back to UTM. Preprocessing refuses force_reproject together with a rotation.
Vector inputs in a frame CRS go in FlatGeobuf or GeoPackage. GeoJSON cannot record a custom CRS, so a GeoJSON written in a frame CRS is read as WGS84. Tools whose help asks for a vector in the DEM's CRS (a crossing line for lsdtt-basin-extraction) need it in the frame CRS; lsdtt-swath-profile reads its baseline's coordinates as-is, which is why the GeoJSON axis above works.
The strip's long sides are grid edges, so boundary conditions can be set per side, for example no flux across the strip and open at its ends.
Nodata wedges need no special handling. Where the source DEM is smaller than the frame, the corners are nodata, and the fills drain cells beside nodata as at any edge.
Lost the frame file? run_rotate_back without rotation_frame_fname rebuilds the frame from each raster's own CRS and rotates back to the UTM zone of the frame centre on the lattice anchored at (0, 0), with a warning: a source DEM origin off that lattice cannot be recovered without the file.

All the parameters are in the lsdtt-fetch-raster and lsdtt-raster-preprocessing references.

How oblique frames work

This section describes how lsdtt3 builds the frames, how rasters are warped into a frame and back, what happens to bearings, and how far results agree between frames and UTM.

The code is in src/core/lsdtt_oblique_frame.{hpp,cpp} (frame geometry and the feature-file reader), src/core/lsdtt_rotation_frame_file.{hpp,cpp} (the frame file), src/core/lsdtt_raster_warp.{hpp,cpp} (warping into a frame and back), src/core/lsdtt_raster_fetch.cpp (fetching into a frame) and src/core/lsdtt_true_north.{hpp,cpp} (grid convergence and true-north bearings).

The projection

A frame is a Hotine Oblique Mercator projection, variant B (PROJ omerc), whose central line runs along the feature:

+proj=omerc +lat_0=<centre lat> +lonc=<centre lon> +alpha=<azimuth> +gamma=90
            +k=1 +x_0=0 +y_0=0 +<datum of the source> +units=m +no_defs

Feature along +x. The axis azimuth is folded into (−90°, 90°] (a line has no direction), and gamma=90 turns the rectified grid so the central line is the x axis. Without the fold, an azimuth of 135° put the line on the y axis. Azimuths within 1e-4° above −90° are taken as +90°: PROJ places the frame centre half a meridian away for alpha within about 1e-6° of −90.
Centre at the origin. The frame centre (the geodesic midpoint of the axis end points) maps to (0, 0), so the extent is symmetric: x in [−L/2 − b, L/2 + b], y in [−W/2, W/2] for axis length L, end buffer b and width W (b also on the long sides for a polygon or a data footprint), snapped outward to whole cells.
Scale. The scale factor is 1 on the central line and grows as about d²/(2R²) at distance d from it: 3 × 10⁻⁵ at 50 km, against 4 × 10⁻⁴ on a UTM zone's central meridian.
Datum. The frame is built on the geographic CRS of its source: WGS84 for a fetch (Copernicus GLO-30 is WGS84), the source DEM's own datum when an existing DEM is rotated, so no datum transformation enters the warp.
Self-check. Every frame is checked when built and when read: its centre must map to (0, 0) within 1 mm, and a point placed along the azimuth (geod_direct through an azimuthal equidistant projection) must land on +x within 0.01 %.

A raster in a frame is an ordinary north-up raster in that CRS, so the DEM loader, the flow routing and every analysis read it unchanged.

Building a frame from a feature

Feature	Axis	Width
Line or multi-part line	the two furthest vertices of all parts	given (`*_width`)
Two points	the two points	given
Polygon	longest side of the minimum-area enclosing rectangle	the rectangle's short side
Valid data of a raster	as a polygon, from the outer corners of the first and last valid cell of every row	the rectangle's short side

The end buffer defaults to ten cells. A polygon whose rectangle has long/short below 1.1 is warned as having an ill-defined axis. Geodesics (lengths, azimuths, midpoints) are ellipsoidal on the frame's datum.

Warping into a frame and back

Into a frame. Fetching warps the source tiles straight into the frame; a DEM on disk is warped from its own CRS. Either way the DEM is resampled once. The cell size is grid_spacing if given; otherwise a DEM projected in metres keeps its own, and any other DEM gets the size reprojection would choose.

Rotate-back target. The frame file records the grid that products go back onto. After rotating a DEM projected in metres, it is that DEM's own CRS, cell size and lower-left corner, so products land on its pixels (a whole number of cells away). After a fetch, or for a geographic source, it is the UTM (or NSIDC polar stereographic beyond 75°) CRS of the frame centre on the 0, 0 lattice that plain fetches use. A raster at the frame's cell size goes to the target cell size; any other keeps its own.

Fixed resampling kernel. GDAL sizes its resampling kernel per warp chunk from the ratio of the chunk's source bounding box to its destination size. A rotated grid's source box is larger than the chunk, so an equal-cell rotation is treated as downsampling and the kernel is widened by an amount that depends on the chunk size, the rotation and the output driver. Frame warps therefore fix the kernel scale (GDAL's XSCALE/YSCALE: 1, or source/target pixel size when the frame cell is coarser) and use the exact transform (-et 0). For a 30 m Copernicus DEM rotated by 36° (5.6 million cells), the warp then agrees with a fixed-scale exact gdalwarp to 0 m whatever the warp memory or output format, where the default kernel had differed from it by 0.87 m median and up to 55 m on cliffs.

Round trips. Warping a DEM into a frame and back (cubic both ways) on that DEM changed elevations by a mean of 0.0002 m, 0.36 m RMS, 1.3 m at the 99th percentile and 18 m at most; on a smooth synthetic DEM, by 0.0007 m median and 0.22 m at most. Categorical rasters rotated with near keep their values: 184 of 184 IDs survived a round trip with no new values.

The rotation frame file

<name>_rotation_frame.txt holds one key: value per line, written so the file is also a flat YAML mapping (yaml.safe_load reads it): text values single-quoted, WKT on one line, numbers with 17 significant digits and always a decimal point, booleans as true/false. It records provenance (frame_file_version, lsdtt3_version, created_utc, created_by, rotated_dem_fname), the frame (centre_lon, centre_lat, azimuth_deg, cell_size, frame_x_min … frame_y_max, base_crs_wkt, frame_crs_wkt), how it was defined (frame_source, feature_fname, feature_layer, the axis end points, axis_length_m, width_m, end_buffer_m, buffer_all_sides) and the rotate-back target (target_crs_epsg or target_crs_wkt, target_cell_size, target_origin_x, target_origin_y, source_dem_fname).

It has its own reader rather than the parameter-file parser, which ends a value at the first space and caps a line at 1023 characters. Reading refuses a missing or newer version, duplicate or missing keys, a frame_crs_wkt that does not match the stored centre and azimuth, and a frame that fails the axis self-check; it warns on unknown keys and when the frame rebuilt from the stored definition differs from the stored one.

Without the file, a frame is rebuilt from a raster's own CRS and extent, and the target falls back to the analysis CRS of the frame centre on the 0, 0 lattice, with a warning: an unsnapped source origin cannot be recovered.

Bearings and grid convergence

Every bearing lsdtt3 writes is clockwise from grid up. In a frame with the feature along +x, 90° is along the feature. The grid convergence c(x, y) is the angle of true north measured clockwise from grid up, computed by stepping 1e-4° north through the CRS's own geographic CRS; a grid bearing g is the true bearing g − c. At a frame's centre c = 90° − azimuth; on UTM, |c| is at most about 3°.

With bearings_relative_to_true_north, each output is converted at its own location: a raster cell by cell (c on a lattice every 32 cells, interpolated bilinearly: within 2e-6° of a per-cell conversion on a 45 km × 12 km frame), a point output at its point, a bearing between two points at their midpoint, a basin or patch axis at its centroid, and a whole-raster 2D spectrum with c at the grid centre (written in the CSV header). Angles between two directions are not converted. The D8 diagnostic fraction_on_d8_multiples in the channel-bearing statistics is always taken from grid bearings. Hillshade illumination is never converted.

Measured against pyproj: c to within 1e-3° in UTM north and south, NSIDC polar stereographic and frames; a plane dipping to true azimuth 200° gives a converted polyfit aspect of 200.000° in UTM and in frames at azimuths −60°, 0°, 30° and 89°.

Agreement between frames and UTM

Analytic surfaces were defined in UTM 12N metres at 38.4° N, −110° (600 × 600 cells of 10 m, centred off every lattice node) and analysed on the UTM grid and on frames at azimuths 0°, 15°, 30°, 45°, 60°, 75° and 89.9°:

Quantity	UTM	Frames, all azimuths
Plane slope (true 0.2), median	0.2	0.199939
Paraboloid Laplacian (true −1e-3 /m), median	−1.000e-3	−0.99939e-3
Funnel D8 outlet area (km²)	19.6357	19.6464 to 19.6472
Funnel D-infinity and MFD outlet area (km²)	19.4361	19.4468 to 19.4476
Funnel D8 path length from the rim / radius	1.0548	1.0550 to 1.0558
V-valley D8 channel length / straight length	1.0008	1.005 (0°, 45°, 89.9°); 1.073 to 1.076 (15°, 30°, 60°, 75°)

Frames agree with each other to better than 4 × 10⁻⁵ in slope, curvature and area. The offset between frames and UTM (−3.1 × 10⁻⁴ in slope, +6.1 × 10⁻⁴ in area) is the UTM scale factor 1° from the zone's central meridian (k = 0.99969), not the rotation. D8 channel lengths are not invariant: they depend on the channel's angle to the grid (+0.5 % along an axis or diagonal, +7.3 to 7.6 % at 15° off), the known D8 path bias, which is the same on UTM for a channel at that angle.

Storage

For a 100 km × 10 km strip at 45° on the Colorado Plateau (Copernicus 30 m), the frame holds 3354 × 334 = 1,120,236 cells (GeoTIFF 3.7 MB); the same strip on a UTM grid needs 2629 × 2604 = 6,845,916 cells (22.0 MB), 6.1 times as many. Channel extraction on the frame peaked at 143 MB resident memory, against 618 MB on the UTM grid.

Figures

lsdviztools3 draws a frame raster in the Cartesian X projection in frame km, so the strip is a rectangle. GMT's own oblique Mercator (-JO) is not used: it is spherical. mapproject -Joa gives the same coordinates whatever ellipsoid GMT is told to use, and matches PROJ's omerc on a sphere to the centimetre; against the ellipsoidal frames lsdtt3 writes, a point 50 km along a 30° axis at 38.4° N is 33 m long and 89 m off the axis (177 m at 100 km). The north arrow is turned by the grid convergence at its own position, and the graticule is the reference CRS's easting/northing grid transformed into the frame.