Fixed-width curvilinear swath profile

At the beginning

Fixed-width curvilinear swath profile is equivalent to traditional swath profile that typically stacks a set of profile lines perpendicular to a baseline within a swath shape, or band, from which the elevation or other types of data (e.g., precipitation, relief, roughness, etc.) are sampled. Runing fixed-width swath analysis in PyOSP follows simple workflow. Here, we use a synthetic case to illustrate the implementation steps.

Figure below shows an artificial mountain range with a homogeneous width (90m) situated above the flat ground (0m).

As suggested by its name, the fixed-width swath analysis requires a defined width which is set as 100 m here in order to encompass the cross width of the mountain (90 m). Additionally, the baseline is set across the longitudinal axis of the mountain. The following shows the code snippet to generate a swath object.

Generate a swath object

import pyosp

baseline = pyosp.datasets.get_path("homo_baseline.shp") # the path to baseline shapefile
raster = pyosp.datasets.get_path("homo_mount.tif")  # the path to raster file

orig = pyosp.Orig_curv(baseline, raster, width=100,
                       line_stepsize=3, cross_stepsize=None)

Processing: [#########################] 71 of 71 lineSteps

Swath polygon and polylines

Using methods out_polygon() or out_polylines(), we can generate shapely objects to be plotted with Matplotlib

PyOSP also provides function read_shape to let us read the baseline shapefile into shapely object, and plot with polygon and polylines as below:

import matplotlib.pyplot as plt
from pyosp import read_shape

# read the baseline shape
line_shape = read_shape(baseline)
lx, ly = line_shape.xy

# Plot the swath polygon
fig, ax = plt.subplots()
swath_polygon = orig.out_polygon()
px, py = swath_polygon.exterior.xy
ax.plot(px, py)
ax.plot(lx, ly, color='C3', label="Baseline")
ax.set_aspect('equal', adjustable='box')
ax.set_title("Swath polygon")
ax.legend()

# Plot the swath profile lines
fig, ax = plt.subplots()
swath_polylines = orig.out_polylines()
for line in swath_polylines:
    x, y = line.xy
    ax.plot(x, y, color='C1')

ax.plot(lx, ly, color='C3', label="Baseline")
ax.set_aspect('equal', adjustable='box')
ax.set_title("Swath profile lines")
ax.legend()

<matplotlib.legend.Legend at 0x7fa4aa7d8700>

Noticing that the polygon is generated by connecting ending points of polylines. In cases that polylines are heavily overlapping with each other, the presented swath polygon could be less reliable. It is worth noting that data is always sampled through the locations of polylines, rather than polygon.

The above also shows that there are cutting conors around the swath area. Now, let’s investigate this by plotting together with raster dataset.

import numpy as np
import gdal

ds = gdal.Open(raster)
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()

xRes = gt[1]
yRes = gt[5]

# get the edge coordinates and add half the resolution
# to go to center coordinates
xmin = gt[0] + xRes * 0.5
xmax = gt[0] + (xRes * ds.RasterXSize) - xRes * 0.5
ymin = gt[3] + (yRes * ds.RasterYSize) + yRes * 0.5
ymax = gt[3] - yRes * 0.5

ds = None

# create a grid of xy coordinates
xy = np.mgrid[xmin:xmax+xRes:xRes, ymax+yRes:ymin:yRes]

fig, ax = plt.subplots()
ax.plot(px, py, color="gold", label='swath polygon')
im = ax.pcolormesh(xy[0], xy[1], data.T, cmap=plt.cm.terrain)
ax.plot(lx, ly, color='C3', label="Baseline")
ax.set_aspect('equal', adjustable='box')
ax.set_xlim([-10, 210])
ax.set_ylim([-10, 210])
plt.colorbar(im)

<matplotlib.colorbar.Colorbar at 0x7fa4aa5463a0>

As shown above, the swath area was bounded by the range of raster dataset. During the calculation of swath object, PyOSP keeps detecting the boundary of dataset to prevent sampling into no data region.

PyOSP can also export generated polygon or polylines as shapefiles. Hence the result can be read and processed in other GIS softwares like ArcGIS, QGIS, etc. For example:

# This will generate polygon and polyline shapefiles in the current folder
pyosp.write_polygon(orig.out_polygon(), "./orig_polygon.shp")
pyosp.write_polylines(orig.out_polylines(), "./orig_polylines.shp")

Swath profile

To plot the corresponding swath profile, simply run method:

orig.profile_plot()

If the figure needs to be customized, an axis argument can be passed to the object:

fig, ax = plt.subplots()
orig.profile_plot(ax=ax, color="maroon")
ax.set_xlabel("Length (m)")
ax.set_ylabel("Depth (m)")
ax.grid()

Summary

This section covers the basic usage of PyOSP to perform traditional fixed width swath profile. For more examples and advance topics, please refer to our example gallery.