Projected surfaces

Most GIS packages report planar area, not surface area. Some packages, such as ArcGIS (3D Analyst) and Surfer provide the option to compute surface area. Computation of surface areas will result in a value that is greater than or equal to the equivalent planar projection of the surface. The ratio between the two values provides a crude index of surface roughness.

If a surface is represented in TIN form, the surface area is simply the sum of the areas of the triangular elements that make up the TIN. Let Tj={xij,yij,zij} i=1,2,3 be the coordinates of the corner points of the jth TIN element. Then the surface area of the planar triangular (or planimetric area) determined when z=constant, can be calculated from a series of simple cross-products:

Note that this expression is the same as the polygon area formula we provided earlier, but for a (n‑1)=3‑vertex polygon:

where i+1=4 is defined as the first point, i=1. The term inside the brackets is the entire area of a parallelogram, and half this gives us our triangle area. This expression can also be written in a convenient matrix determinant form (ignoring the j’s for now) as:

If the slope of Tj, θj, is already known or easily computed, the area of Tj can be directly estimated as Aj/cos(θj) (Figure 4‑4). Alternatively the three dimensional equivalent of the matrix expression above can be written in a similar form, as:

If z=0 for all i, the second and third terms in this expression disappear and we are left with the result for the plane. Suppose that we have a plane right angled triangle with sides 1,1 and √2, defined by the 3D coordinates: (0,0,0), (1,0,0) and (1,1,0). This plane triangle has area of one half of a square of side length 1, so its area is 0.5 units, as can be confirmed using the first formula above. If we now set one corner, say the point (1,1,0) to have a z-value (height) of 1, its coordinate becomes (1,1,1) and the area increases to 0.7071 units, i.e. by around 40%.

Figure 4‑4 Planimetric and surface area of a 3D triangle

If the surface representation is not in TIN form, but represented as a grid or DEM, surface area can again be computed by TIN-like computations. In this case the 8-position immediate neighborhood of each cell can be remapped as a set of 8 triangles connecting the center-points of each cell to the target cell. Connecting all cell centers in this way creates a form of TIN and surface computation proceeds as described above – this is the procedure used in ArcGIS (3D Analyst). A more sophisticated method involves a similar form of computation but calculated on a cell-by-cell weighted basis. In this procedure the surface area for each cell is taken as a weighted sum of eight small triangular areas. The area of triangles centered on the target cell are computed and 25% of each such area is assigned to the central cell (shaded gray in Figure 4‑5B).

Figure 4‑5 DEM surface area: A. 3x3 DEM element

Figure 4‑5 DEM surface area: B. Plan view of triangular components

Research has suggested that the resulting total surface area for the cell is close to that produced by TIN models for grids of 250+ cells. We can examine this model more closely by taking a specific example. Let us assume that the 3x3 DEM surface segment in Figure 4‑5 has cells which are 1 meter square and elevation values given in meters as shown below:

The nine individual cells can be referenced in various ways – below use a simple numerical referencing around the central or target cell:

The planar area of each cell is simply 1x1, so 1 square meter. To compute the estimated surface area we need to calculate the area of each of the 8 sloping triangles and allocate 25% of this area to the central cell. To do this we use the determinant expression shown earlier. For the triangle shown in Figure 4‑5A (the triangle connecting cells 1-2-0) this gives an area of 3/2, so 25% of this gives 3/8 units as opposed to the 1/8 units that to a flat triangle such as 0-4-5. Continuing this computation for all 8 triangles gives the overall estimated surface area using this algorithm as 2.5 sq.m. The process is then repeated for every pixel, adjusting for edge effects where necessary as with many other operations of this type (e.g. by adding a temporary row or column to each edge that duplicates the existing values).

This example is unrealistic of course since In most DEM datasets the horizontal and vertical scales are not the same. For example, if the individual cell sizes in the example above are 10 meters x 10 meters and the vertical measurements are in meters, the procedure yields a surface area approximately 3% greater than the planar area. Grid surface areas can also be computed using the slope adjustment method described above assuming that a slope grid has been first been computed, or using grid data with linear interpolation to facilitate simpler 3D quadrilateral computations.

Surface areas are computed relative to a reference plane, usually z=zmin or z=0. In some instances it may be necessary to specify the value to be associated with the reference plane, and whether the absolute values, positive values and/or negative values are to be computed. Surfer, for example, supports all three options which it describes as positive planar and surface area, negative planar and surface area and total planar and surface area. For example, consider the surface illustrated in Figure 4‑6. This shows a surface model of a GB Ordnance Survey DEM for tile TQ81NE, which we shall use as a test surface in various sections of this Guide. This is a 5000mx5000m area provided as a 10mx10m grid of elevations, with an elevation range from around 10m‑70m. Using a reference plane of z=30m the surface area above this reference plane is roughly 1.3sq kms and below 30m is roughly 2.1sq kms (3.4sq kms in total compared to 2.5sq kms for the planar area). For details of volume computations, which also use reference surfaces, see Section 6.2.6, Pit filling.

Figure 4‑6 Surface model of DEM

Terrestrial (unprojected) surface area

For large regions (e.g. “rectangular” regions with sides greater than several hundred kilometers) surface areas will be noticeably affected by the curvature of the Earth. Such computations apply at large State and continental levels, or where the areas of ocean, ice or cloud cover are required over large regions. Computation of areas using spherical trigonometry or numerical integration is preferable in such cases.

The area of a spherical rectangular region defined by fixed latitude and longitude coordinates is greater than the area bounded by lines of constant latitude — the latter are not great circles, although the lines of longitude are. This effect increases rapidly nearer the poles and is weakest close to the equator. The MATLab Mapping Toolbox (MMT) function AREAQUAD provides a convenient way of computing the area of a latitude-longitude “quadrangle”, although it can be computed without too much difficulty using spherical trigonometry. The area north of a line of latitude, φ, is the key factor in the calculation (the surface area of an entire sphere is 4πR2). The former area is:

where R is the radius of the spherical Earth, e.g. 6378.137 kms. The area of the quadrangle is then simply the difference between the areas, A1 and A2, north of each of the two lines of latitude, φ1 and φ2, adjusted by the proportion of the Earth included in the difference in longitude values, λ1 and λ2:

A (simple) polygon with n vertices on the surface of the Earth has sides that are great circles and area:

where again R is the radius of the Earth and the θi are the internal angles of the polygon in radians measured at each polygon vertex. The simplest spherical polygon is a triangle (n=3 so n‑2=1) and on a sphere the sum of its internal angles must be greater than 180 degrees (π), so the formula equates to 0 when the triangle is vanishingly small. A spherical triangle with internal angles that are all 180 degrees is the largest possible, and is itself a great circle. In this case the formula yields A=2πR2, as expected.

All polygon arcs must be represented by great circles and internal angles are to be measured with respect to these curved arcs. The term in brackets is sometimes known as the spherical excess and denoted with the letter E, thus A=ER2. Note that every simple polygon on a sphere divides it into two sections, one smaller than the other unless all vertices lie on a great circle. If only the latitude and longitude values are known, the internal angles must be computed by calculating the true great circle bearing (or azimuth) from each vertex to the next — this can be achieved using spherical trigonometry or, for example, the MMT function AZIMUTH, which supports explicit geoid values. Azimuth values are reported from due north, so computations of internal angles need to adjust for this factor. Alternatively the MMT function AREAINT may be used, ideally with high point density (i.e. introducing additional vertices along the polygon arcs). This latter function computes the enclosed area numerically. The open source GIS, GRASS, provides a number of similar functions, including the source code modules area_poly1, which computes polygon areas on an ellipsoid where connecting lines are grid lines rather than great circles, and area_poly2 which computes planimetric polygon area.