The 3D Hydrography Program (3DHP) is the first systematic remapping of the hydrography of the US since the original USGS 1:24,000-scale topographic mapping program. The core goal of the program is the production of a seamless, high-resolution hydrographic dataset that is derived from and spatially and temporally integrated with the elevation data collected by the 3D Elevation Program (3DEP).

When elevation-derived hydrography (EDH) data for a drainage basin is submitted for evaluation as candidate 3DHP data, we perform a suite of validation checks, many of which assess how well the data is integrated with the source digital elevation model (DEM) by performing terrain analysis.

We use implementations of terrain analysis algorithms from open source software packages such as Whitebox Tools and SAGA, and we use the Python open source geospatial stack including rasterio, NumPy, and SciPy to create our own implementations. A limitation that these implementations share is that they tend to work with the entire dataset in memory, requiring multiples of the DEM's uncompressed size in RAM to complete. For elevation datasets that may exceed 100 GB, provisioning these systems would be impractical and increasingly cost-prohibitive.

The standard approach to this problem is tiling – splitting the large dataset into overlapping rectangular chunks that can be processed individually. This works well with local functions or neighborhood functions with a fixed or maximum neighborhood size, in which cases it is possible to eliminate edge effects by using a sufficiently large tile overlap. In the case of flow accumulation algorithms, however, edge effects are unavoidable, as values may accumulate from one end of an elevation dataset to the other. For zonal statistics, too, features may extend beyond any reasonable tile overlap limit.

In this presentation we will discuss the solutions we use for these edge-case problems of raster processing with larger datasets. We will define our unique requirements, provide an overview of the technical details of the solutions, and discuss the benefits and limitations of these approaches and how they may or may not be applicable to other situations.

To minimize edge effects in flow accumulations, we developed a watershed-based tile generation process that uses GeoPandas, shapely, and GDAL. Water accumulates throughout a watershed and flows to the next watershed downstream at a pour point, so by locating tile boundaries at watershed boundaries, the edge effects in flow accumulation patterns are only pronounced at watershed inflows and outflows. Including tile overlap tends to minimize these effects because channels are generally well-defined at the pour point of a watershed. We used the 12-digit Watershed Boundary Dataset (WBD) boundary polygons as the initial tile boundaries, as they are the smallest watershed boundaries available nationally. The tile generation process computes the anticipated raster tile dataset size for each watershed from the DEM resolution and bit depth, splits the watershed if necessary to achieve the target tile size, then merges adjacent tiles iteratively to minimize the total number of tiles and area of overlap. We will discuss the details of this algorithm in the presentation.

For zonal statistics of vector features and other raster manipulations that relate to vector features, we leveraged rasterio's geometry windows to read and write elevation data for individual features, which not only minimized memory usage but also enabled parallel execution and fine-grained control of analysis steps. In the presentation we will discuss how to implement this pattern and how we applied it to perform zonal statistics, hydro-enforcement, and elevation profiling. We will also discuss the limitations of this approach, the major one being that features in some datasets may be large enough to cause memory issues.