Data preparation involved several manual steps. This included downloading the FWA spatial files, merging the stream network feature classes into a single layer, and converting to a GeoPackage layer. The FWA river polygons and coastline were also converted to layers within the same GeoPackage. River polygons were further prepared by running the Dissolve Boundaries tool in ArcPro. This decision was made so that each river polygon could be assigned the maximum stream order value from the FWA lines within it. Without dissolving, when a river splits near the mouth, one or more splits will have a very low order while one split will retain a higher order. By dissolving the shared boundaries and assigning the maximum order value to that polygon, this problem is addressed. The Fraser river dissolved polygon was buffered by 6 km and manually edited to extend into the delta and generally align with modelled salinity data. Existing simplified SDM boundary feature classes were also copied into the GeoPackage.

The following steps describe the geoprocessing operations that are part of the code, which was run separately for each of the marine regions, with the same arguments provided to the script. A bounding box was created in memory by loading a bathymetry raster from the SDM predictor layers. The FWA coastline was also loaded into memory. The bounding box was used as a spatial filter when loading other vector datasets, for space and time efficiency. A regional simplified SDM boundary vector dataset was loaded into memory and used to clip the FWA stream network and dissolved river polygons. Geometry was checked and converted to 2D if it had Z coordinates.

Points were created by intersecting the clipped FWA polylines and FWA coast polylines (retaining only the points of intersection). Any duplicate points (based on geometry) were removed. These points represent river mouths.

A spatial join operation was run on the river mouth points and river polygons in order to attach the required stream variable (i.e., STREAM_ORDER) to the rivers. After the join, river polygons were clipped by the simplified SDM boundary for time efficiency. Points were then created within each of the river polygons. This was done by iterating the river features, creating a rectangular mesh of points withing the bounds of the clipped river geometry, and removing points not contained by the actual geometry of the river polygon. These points were added to the existing river mouth points with the stream variable values. Additional points were added to the GeoDataFrame along line segments of the FWA stream network clips at ½ the cell size of the bathymetry layer.

The points are written to disk as a layer in an output GeoPackage. Then, the points in memory are converted to a raster, using the integer STREAM_ORDER values for the values in the raster (NoData is represented by 0 for the output raster). This single-band raster is written to disk as a GeoTIFF.

The stream variable raster was read into memory as a numpy array and the values were transformed by applying a power function (raising the values by a factor of 0.5) in order to decrease the range of values. A focal window (np.meshgrid) of weights was created for each bin (number of unique STREAM_ORDER values) where the cells at the centre have the greatest value, and the cells at the edges are 0. For each stream order, a different sized window was used. For example, stream order 7 used a larger window than stream order 1. To ensure that the window size for a given stream order was the same between regions, 10 bins were created (Fraser river has stream order 10 which was the max of any regions) and scaled using numpy's geomspace function that return numbers spaced evenly on a log scale (a geometric progression). The stream order window sizes are as follows: {1: 51, 2: 63, 3: 77, 4: 95, 5: 117, 6: 143, 7: 175, 8: 215, 9: 265, 10: 325}. For each stream order, the array was subset by that particular order value, then convolved using ndi.convolve (https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.convolve.html) over the existing array with a unique window size. On each iteration, the existing convolved array was summed with the current convolved array (additive process during each iteration). Once complete, the masked and convolved array was written to disk as a GeoTiff.

After all regions were completed individually, a separate script (https://gitlab.com/dfo-msea/environmental-layers/normalize-fetch) was run to rescale (normalize) the layers against each other. In this way, the cell values from any given region represent the same relative influence between regions. If viewing the data in GIS software, setting the stretch symbology type to ‘Histogram Equalize’ will display the more local changes more obviously than some of the other symbology types and may be useful where the range in values is greater (i.e. the Salish Sea region).