The recent surge of applications using terrestrial cosmogenic nuclides (TCNs) to calculate catchment-averaged erosion rates from isotopic concentrations in fluvial sediment, and the prospect of coupling TCN production functions with numerical surface process models (SPMs), necessitate a fast and accurate algorithm for the calculation of topographic shielding. Topographic shielding refers to the proportion of the incoming cosmic radiation that is shielded by the surrounding topography, the scaling factor being defined as the ratio of the unshielded (total minus shielded) to the total (or maximum) cosmic ray flux (i.e. the flux received by a horizontal, unobstructed surface). Topography contributes to the reduction of TCN production by obstructing a certain proportion of the incoming flux and by modifying the angle of incidence. Available algorithms calculate the proportion of obstructed radiation by dividing the horizon as seen by the sample (a grid cell in the case of a DEM), into arc segments (usually of equal length) for which the average obstruction heights expressed as zenith angles are calculated. The use of these methods is feasible only when dealing with a small number of isolated samples, since the identification of obstructions when dealing with an entire area is computationally very intensive. This paper describes a method that uses a relief shadow modelling technique to identify those areas of a DEM that are under shadow (i.e. shielded), and thus to account for the obstructed radiation. This method produces results that are very similar to those obtained using a direct implementation of available methods (maximum difference between results of c. 0·1). The method based on relief shadow modelling is also faster than a direct implementation of any available method and can be readily implemented in any GIS system with raster capabilities.