spatial analyst - Sum raster data visible from given point, by proportion of view

I am wondering if there is a tool that does the following (preferably in Arc but other options considered).

Simple example: what proportion of the view from a given point on the landscape is occupied by a given land use type?

You might think you could do this by computing a viewshed and using zonal statistics to sum the land use over the area of the viewshed. However, this will give an answer as a proportion of land area, which will therefore be biased towards more distant land uses - if you can see a long way, you can see an enormous area even if it only occupies a small proportion of your view. So as far as I can tell, the tools in Arc do not do the job.

More complex example: Sum the weights on a raster by proportion of view they occupy from a given point. Essentially the same as the first example but now dealing with continuous values rather than 1 and 0 identifying the land use type.

More complex example yet: Sum the weights on a raster by proportion of view they occupy for each point in the output raster. Same as the first example but repeated for every point.

Answer

The inaccuracies introduced by ignoring buildings, trees, etc., suggest that we don't need extremely high accuracy in this computation. If we intend to exclude the cell on which the observer rests, along with immediately neighboring cells, then the calculations can be greatly simplified. This answer describes the simplified calculations.

The view subtended by any given cell is proportional to all of the following:

• 
  

  The cell's actual area in three dimensions. This accounts for the slope of the land over that cell; steeper slopes indicate greater areas in proportion to the secant of the slope.

  
  


• 
  
  

  The cell's distance. Because a cell cannot have a constant distance from an observer, this is an approximation. For cells reasonably distant, we may use the distance between the observer and the cell's center to approximate distances throughout the cell.

  
  


• 
  

  The sine of the angle made between the cell and the observer. More specifically, the inner product of the normal direction to the cell and the normalized vector towards the observer gives the relative proportion of the cell's area subtended in the observer's field of view.

  
  

To compute this inner product, we need to obtain the x, y, and z coordinates of two vectors: (1) a vector pointing from the cell's center back toward the observer and (2) a vector perpendicular to the surface at the cell itself. (1) is perhaps most easily obtained from two components, the horizontal and vertical. The vertical component is the difference between the cell's elevation and the observer's (total) elevation. The horizontal component equals the product of the distance from the cell to the observer and the direction vector (that is, cosine and sine) of the cell towards the observer. Why decompose it like this? Because the angle from the cell towards the observer is given by the aspect of the Euclidean distance grid relative to the observer. (This is obvious once you recognize that the Euclidean distance describes an abstract "surface" that is an inverted code centered at the observer.)

We obtain the (x,y,z) coordinates for (2) from the slope and aspect of the DEM at the cell itself.

Here are the details, with a running example.

1. 
   

   We begin with a DEM whose values are in the same units of measure as its horizontal coordinates. The white dot shows the location of an observer whose height is slightly above the DEM itself.

   
   

   

   
   


2. 
   

   Compute the slope (in degrees or radians, not percent) of the DEM and its aspect.

   
   


3. 
   

   Compute the Euclidean distance grid relative to the observer.

   
   



4. 
   

   Compute the aspect of the distance grid (in degrees or radians).

   
   


5. 
   

   Compute the (x,y,z) coordinates of the vectors pointing from the cell centers back toward the observer. In peudocode, assuming all angles are in degrees, these computations are

   
   

   X = Cos((90.AsGrid- [Aspect of Distance]) / 180 * 3.14159265358979323) * [distance]
   Y = Sin((90.AsGrid- [Aspect of Distance]) / 180 * 3.14159265358979323) * [distance]
   Z = (Observer's total elevation) - [DEM]
   

   
   

   For the latter, you need to obtain the observer's elevation by interpolating the DEM's height at the observer's coordinate and adding the observer's height above the DEM. The sum is a number to be used for Observer's total elevation.

   
   
   


6. 
   

   The (x,y,z) coordinates of the unit normal vector at each cell are obtained by the (usual) formulas converting between spherical and 3D Cartesian coordinates: 90 Degrees - [Aspect] is spherical longitude and 90 Degrees - [Slope] is spherical latitude. Whence

   
   

   Nx = Sin([Aspect] / 180 * 3.14159265358979324) * Sin([slope] / 180 * 3.14159265358979324)
   Ny = Cos([Aspect] / 180 * 3.14159265358979324) * Sin([slope] / 180 * 3.14159265358979324)
   Nz = Cos([slope] / 180 * 3.14159265358979324)
   

   
   

   Here, for instance, is the [Nx] grid, with red cells positive and gray cells negative:

   
   

   

   
   



7. 
   

   Compute the relative area of each cell as the reciprocal of the cosine of the slopes.

   
   


8. 
   

   The distance in three dimensions between each cell's center and the observer is given by the Pythagorean theorem as

   
   

   [Distance 3D] =  (([X] * [X]) + ([Y] * [Y]) + ([Z] * [Z])) ^ (1/2)
   



9. 
   

   Compute the weights. As described, this is (a) the cell's area times (b) the inner product of the normal and distance vectors, (c) normalized by the length of the distance vector, and (d) divided by the square of the distance. The net effect of (c) and (d) is to divide by the cube of the distance:

   
   

   [Weight] = [Area] * (( [Nx] * [X]) + ([Ny] * [Y]) + ([Nz] * [Z])) / ([Distance 3D]^3).
   
   



10. 
    

    Exclude negative weights using SetNull or Con: negative weights correspond to cells that will not be visible because they slope away from the line of sight from the observer. Here, the result is shown on a logarithmic scale (yellow=large, blue=small, black=invisible):

    
    

    

    
    

This grid should look exactly as if a light were shining from the observer's location onto the DEM's surface: the values on the grid represent the intensity of that light at each cell, which is proportional to the view subtended by the cell from the observer's point of view.

With these weights, every version of this problem can be directly solved. Given a raster of any kind of values, and given another raster indicating which values should be used (such as an indicator of visibility), it is desired to form the weighted sum of the values in the first raster. (It is a good idea to exclude very nearby cells from this calculation, due to the inaccuracies of the approximations there.) This is obtained with two zonal sums:

1. 
   
   

   Compute the zonal sum of the weights, masked by the visibility grid.

   
   


2. 
   

   Compute the zonal sum of the product of the weight and value grids, masked by the visibility grid.

   
   

Dividing the result of (2) by the result of (1) yields the weighted average of values within the mask.

This procedure needs to be repeated separately for each observer (by looping over the observers). There's no way around this: every set of weights is observer-specific.