| GISdevelopment.net ---> Technology ---> GIS |
Surface area processing in GIS
Sandeep N Kundu
IFGI, University of Münster, 48149 Münster, Germany
Phone: 0049 251 8330058 Fax: 0049 251 8339763
Email: skundu@ifgi.uni-muenster.de
Biswajeet Pradhan
Asian Institute of Medicine, Science & Technology,
08000, Sungai Petani, Kedah Darul Aman, MALAYSIA
Phone: 00604 445 2233 Ext. (222) Fax: 00604 442 2887
Email: biswajeet@aimst.edu.my (presenting author)
Sandeep N Kundu
IFGI, University of Münster, 48149 Münster, Germany
Phone: 0049 251 8330058 Fax: 0049 251 8339763
Email: skundu@ifgi.uni-muenster.de
Biswajeet Pradhan
Asian Institute of Medicine, Science & Technology,
08000, Sungai Petani, Kedah Darul Aman, MALAYSIA
Phone: 00604 445 2233 Ext. (222) Fax: 00604 442 2887
Email: biswajeet@aimst.edu.my (presenting author)
Abstract
Geographical Information Systems (GIS) provides a framework for modelling real world at different map scales through projection systems and spatial data models (raster or vector). Processing functionalities have been deemed fundamental to such a modelling environment and are further being developed for applications on a real-time basis (Berry, 1987). Distance and area processing in GIS play an important role in the decision making process. It is therefore important to evaluate techniques by the GIS model, to obtain reliable distance and area measurements. Traditionally, distance and area measurements from field surveying methods were entered in the feature attribute table. It is desirable to replace this time and cost consuming practice through accurate processing of spatial features in GIS. The current article discusses and compares techniques for real world area approximation within GIS.
Introduction
Area is a fundamental parameter derived from terrain analysis which is necessary for numerous decision-making processes. There is a need to obtain an accurate approximation of distance and area which are the fundamental terrain parameters associated with spatial data content since these have a bearing on cost (Navigation, cadastral applications). This problem has been addressed within the framework of a given data scale and a suitable projection. In most elementary GIS models, distance and area calculations are based on the vector data model which is a planar approximation of data layers. This does not approximate distance or area on an undulating surface as in reality. Given the pre-conditions that the data layers are in the desired scale and projection, a novel way to integrate a raster elevation model to retrieve the near true area and distance (involving elevation data & slopes) is proposed.
Methodology
The technique utilises a Digital Elevation Model (DEM) layer in addition to the vector polygon layers and experiments on some algorithms for surface area calculations for the polygon parcels. The technique can also be extended to distance calculations.
A DEM is essentially, “a regular gridded matrix representation of the continuous variation of relief over space” (Burrough & McDonnell, 1987). This in fact means that a DEM is a continuous grid containing the elevation at its spatial location in each grid cell. Before the 80s, topographic contour maps and aerial photographs were the primary sources of geomorphological information. The accuracy of these maps was highly questionable and there were also discrepancies in information content for two maps for the same region but from different sources. After the advent of DEM this discrepancy was reduced to a large extent. Nevertheless, a DEM is also subject to error depending of the method of its acquisition as each method has its own advantages and limitations. Usually a DEM is interpolated from digitising spot heights as recorded on field or from an existing contour map. But a better DEM is extracted from stereo aerial photographs or stereo digital remote sensing imagery. A DEM was preferred to a TIN because the latter is based on slope inflection points and fail to capture the subtle difference in slope along a particular aspect.
Calculating Slope from DEM
Onorati, et al. (1992) evaluated the efficiency of several slope calculation methods. The methods they took into consideration were
The slope function in Arc/Info follows the second derivative method. Conceptually, the slope function fits a horizontal plane to the z (elevation) values of a 3 by 3 cell neighbourhood around the processing or centre cell. The direction the horizontal plane faces is the aspect for the processing cell. The slope for the cell is calculated from the 3 by 3 neighbourhood using the average maximum technique (Burrough & McDonnell, 1998). If there is a cell location in the neighbourhood with a no data z value the z value of the centre cell will be assigned to the location. At the edge of the grid, at least three cells (outside the grid's extent) will contain no data as their z values. These cells will be assigned the centre cell's z value. The result is a flattening of the 3 by 3 horizontal that is fit to these edge cells, which thus usually leads to a reduction in the slope. This flattening effect did not affect the current study, as the edge cells were not included in the area of interest. For simplicity and conventional reasons we adopt the Arc/Info methodology for slope.
Despite the wide use of GIS for several terrain analyses, including hydrological and agricultural studies, a basic standardised functionality is still lacking in modern commercial GIS packages, which allows the surface area of a terrain to be calculated from a DEM. Slope gradient being an approximation from the neighbouring pixels, its accuracy in determining the surface area but only remains to be verified. On the other hand, since it is near impossible to determine accurately the actual surface area of a region both because of the fact that traditional methods are not feasible Slope Angle, on a smaller scale and that the surface of a terrain is subject to continuous undulations as compared to the featureless plane. This is the basic reason why the performance of using slope for surface area calculation cannot be intrinsically evaluated. Only a comparison of its performance with respect to other algorithms can be performed up to certain level of satisfaction.
Surface Area
Surface area calculations have been attempted before (Strahler, 1952) and several algorithms were developed for the purpose, all of which are based on slope. Surface area, therefore is a second derivative of elevation data. Elghazali et al. (1986), described the areal parameter, which essentially is a global function that produces a ratio between the surface area and plan area. This parameter was used for terrain characterisation but we found it useful as an measure to estimate surface area.
Surface Area from Slope Gradient The slope area for each pixel is calculated and then the summation of all pixels falling within a parcel constitutes the surface area of the parcel. The slope area for each pixel can be approximated from the resolution of the pixel and the slope value for the pixel. In figure 1, area ABEF is the slope area for the pixel ABCD. The slope area for this pixel is given by the expression (AB*CD)/Cos(è ).
Figure 1: Slope area for a Pixel ABDC
Surface Area from SHR (surface to horizontal area ratio)
Extending the areal parameter, SHR (Surface to Horizontal area Ratio) was conceived, which has a potential to accurately approximate surface area. The algorithm developed was a 3 by 3 roving window operator that triangulates the pixels in the window into eight triangles and then extracting the elevation values from the nine pixels calculates the surface area for each of the 8 triangular faces (Fig. 2).
Figure 2: Pixel base (left) & Surface (right) in a DEM
The equivalent horizontal area was then calculated for the equivalent region within the window (i.e. 4 times the pixel area) and then its ratio to the former (i.e. the SHR value) was stored in the central pixel.
Figure 3: The 3 by 3 roving window for SHR calculation
Figure 4: The DEM and the watershed basin
This operator was programmed in JDK1.3 and was run on the available DEM of the region. The program operates on a 3 by 3 window basis and calculates the surface area of the window by the triangulation (Fig. 3). The ratio (SHR) actually is stored as a percentage so as to minimise errors in float to integer conversions.
SHR = {(SA)/(HA)}× 100
Where SA is summation of surface area for each 8 triangles (Fig. 2) & HA is 4 times the pixel resolution area
The surface area (SAP) for each pixel was calculated from the SHR by the equation below
SAP = {SHR × Res}/100
Where Res is the area of each pixel for the Raster
The value of SHR can never be less than one since the surface area is greater than or equal to the horizontal area. Moreover the upper limit of this value though theoretically is infinity but practically it is finite since a slope value of 90 0 is not practically derivable from a DEM. SHR, being the focal function equivalent of the areal parameter (which is a global function), is suitable for further overlaying operations (local, focal, zonal or global). This was essence of its conception in this research.
Test Case
The algorithms discussed above were applied on a given zone in the study region and a statistical-based comparative study was undertaken to determine their performance and suitability. Surface area as calculated for a delineated watershed basin for the area (Fig. 4) using the aforesaid discussed algorithms were used for pixel based comparative study.
The test area (for which a DEM was readily available) was chosen for its rough topography depicted by a greater range for slope (0 0 -79 0 ).
Results
Summarising the surface area value for the watershed basin we find the SHR method reflecting a higher value as compared to the slope method (Table 1).
Table 1 Surface Area calculated from Slope and SHR
The maximum area for a pixel using the slope method was much higher compared to the SHR method. This indicated that the SHR method was not sensitive to extreme slope values, which is natural owing to its dependency on neighbouring pixel values. Pixels with 0 slope value should ideally yield an SHR value of 1, which is actually not evident from the first record of Table 2. The reason is again attributed to the fact that SHR for 0 slope pixels were dependent on slope of neighbouring cells.
From our studies we could comprehensively conclude that SHR could be used to calculate surface area more accurately for a region with medium to low gradient. This finding makes it amenable that SHR method is better suited to regional and global scale studies whereas slope method can be applied to detailed scale (local) studies. SHR is important for terrain analysis where slope dependency between adjacent pixels is needed to preserve.
Discussion
GIS processing and decision making can be improved utilizing this surface area calculating functionality, especially for calculation of actual length and surface area for objects. Such functionality does add a new dimension to the accuracy of the surface area value although the dependency of accuracy also depends on the elevation model (TIN, DEM, DTED) in use, its information source (interpolated contours, stereo-imagery) and the resolution. Here we provide a comparison between two algorithms using terrain conditions that are amenable for both use. This comparison can also be extended to the study of intrinsic terrain roughness particularly for discriminating terrain types with multi-parameter terrain roughness estimates.
Slope can also be used to calculate the actual length of road and rail network as well, although a special algorithm can be devised to determine the departure of the object from one pixel to another. This departure (similar to aspect) can be used to calculate the apparent slope along that departure and from this apparent slope the actual length of the road or rail network can be approximated. This would enhance the present performance of area and length calculations in GIS and thereby enable decision making in applications involving cost associated area-length calculations.
References
Geographical Information Systems (GIS) provides a framework for modelling real world at different map scales through projection systems and spatial data models (raster or vector). Processing functionalities have been deemed fundamental to such a modelling environment and are further being developed for applications on a real-time basis (Berry, 1987). Distance and area processing in GIS play an important role in the decision making process. It is therefore important to evaluate techniques by the GIS model, to obtain reliable distance and area measurements. Traditionally, distance and area measurements from field surveying methods were entered in the feature attribute table. It is desirable to replace this time and cost consuming practice through accurate processing of spatial features in GIS. The current article discusses and compares techniques for real world area approximation within GIS.
Introduction
Area is a fundamental parameter derived from terrain analysis which is necessary for numerous decision-making processes. There is a need to obtain an accurate approximation of distance and area which are the fundamental terrain parameters associated with spatial data content since these have a bearing on cost (Navigation, cadastral applications). This problem has been addressed within the framework of a given data scale and a suitable projection. In most elementary GIS models, distance and area calculations are based on the vector data model which is a planar approximation of data layers. This does not approximate distance or area on an undulating surface as in reality. Given the pre-conditions that the data layers are in the desired scale and projection, a novel way to integrate a raster elevation model to retrieve the near true area and distance (involving elevation data & slopes) is proposed.
Methodology
The technique utilises a Digital Elevation Model (DEM) layer in addition to the vector polygon layers and experiments on some algorithms for surface area calculations for the polygon parcels. The technique can also be extended to distance calculations.
A DEM is essentially, “a regular gridded matrix representation of the continuous variation of relief over space” (Burrough & McDonnell, 1987). This in fact means that a DEM is a continuous grid containing the elevation at its spatial location in each grid cell. Before the 80s, topographic contour maps and aerial photographs were the primary sources of geomorphological information. The accuracy of these maps was highly questionable and there were also discrepancies in information content for two maps for the same region but from different sources. After the advent of DEM this discrepancy was reduced to a large extent. Nevertheless, a DEM is also subject to error depending of the method of its acquisition as each method has its own advantages and limitations. Usually a DEM is interpolated from digitising spot heights as recorded on field or from an existing contour map. But a better DEM is extracted from stereo aerial photographs or stereo digital remote sensing imagery. A DEM was preferred to a TIN because the latter is based on slope inflection points and fail to capture the subtle difference in slope along a particular aspect.
Calculating Slope from DEM
Onorati, et al. (1992) evaluated the efficiency of several slope calculation methods. The methods they took into consideration were
- The Four Contiguous Right Triangles (FCRT) method,
- Maximum Downward Gradient (MDG).
- Bicubic Spline First Derivatives (BSFD).
The slope function in Arc/Info follows the second derivative method. Conceptually, the slope function fits a horizontal plane to the z (elevation) values of a 3 by 3 cell neighbourhood around the processing or centre cell. The direction the horizontal plane faces is the aspect for the processing cell. The slope for the cell is calculated from the 3 by 3 neighbourhood using the average maximum technique (Burrough & McDonnell, 1998). If there is a cell location in the neighbourhood with a no data z value the z value of the centre cell will be assigned to the location. At the edge of the grid, at least three cells (outside the grid's extent) will contain no data as their z values. These cells will be assigned the centre cell's z value. The result is a flattening of the 3 by 3 horizontal that is fit to these edge cells, which thus usually leads to a reduction in the slope. This flattening effect did not affect the current study, as the edge cells were not included in the area of interest. For simplicity and conventional reasons we adopt the Arc/Info methodology for slope.
Despite the wide use of GIS for several terrain analyses, including hydrological and agricultural studies, a basic standardised functionality is still lacking in modern commercial GIS packages, which allows the surface area of a terrain to be calculated from a DEM. Slope gradient being an approximation from the neighbouring pixels, its accuracy in determining the surface area but only remains to be verified. On the other hand, since it is near impossible to determine accurately the actual surface area of a region both because of the fact that traditional methods are not feasible Slope Angle, on a smaller scale and that the surface of a terrain is subject to continuous undulations as compared to the featureless plane. This is the basic reason why the performance of using slope for surface area calculation cannot be intrinsically evaluated. Only a comparison of its performance with respect to other algorithms can be performed up to certain level of satisfaction.
Surface Area
Surface area calculations have been attempted before (Strahler, 1952) and several algorithms were developed for the purpose, all of which are based on slope. Surface area, therefore is a second derivative of elevation data. Elghazali et al. (1986), described the areal parameter, which essentially is a global function that produces a ratio between the surface area and plan area. This parameter was used for terrain characterisation but we found it useful as an measure to estimate surface area.
Surface Area from Slope Gradient The slope area for each pixel is calculated and then the summation of all pixels falling within a parcel constitutes the surface area of the parcel. The slope area for each pixel can be approximated from the resolution of the pixel and the slope value for the pixel. In figure 1, area ABEF is the slope area for the pixel ABCD. The slope area for this pixel is given by the expression (AB*CD)/Cos(è ).
Figure 1: Slope area for a Pixel ABDC
Surface Area from SHR (surface to horizontal area ratio)
Extending the areal parameter, SHR (Surface to Horizontal area Ratio) was conceived, which has a potential to accurately approximate surface area. The algorithm developed was a 3 by 3 roving window operator that triangulates the pixels in the window into eight triangles and then extracting the elevation values from the nine pixels calculates the surface area for each of the 8 triangular faces (Fig. 2).
Figure 2: Pixel base (left) & Surface (right) in a DEM
The equivalent horizontal area was then calculated for the equivalent region within the window (i.e. 4 times the pixel area) and then its ratio to the former (i.e. the SHR value) was stored in the central pixel.
Figure 3: The 3 by 3 roving window for SHR calculation
Figure 4: The DEM and the watershed basin
This operator was programmed in JDK1.3 and was run on the available DEM of the region. The program operates on a 3 by 3 window basis and calculates the surface area of the window by the triangulation (Fig. 3). The ratio (SHR) actually is stored as a percentage so as to minimise errors in float to integer conversions.
SHR = {(SA)/(HA)}× 100
Where SA is summation of surface area for each 8 triangles (Fig. 2) & HA is 4 times the pixel resolution area
The surface area (SAP) for each pixel was calculated from the SHR by the equation below
SAP = {SHR × Res}/100
Where Res is the area of each pixel for the Raster
The value of SHR can never be less than one since the surface area is greater than or equal to the horizontal area. Moreover the upper limit of this value though theoretically is infinity but practically it is finite since a slope value of 90 0 is not practically derivable from a DEM. SHR, being the focal function equivalent of the areal parameter (which is a global function), is suitable for further overlaying operations (local, focal, zonal or global). This was essence of its conception in this research.
Test Case
The algorithms discussed above were applied on a given zone in the study region and a statistical-based comparative study was undertaken to determine their performance and suitability. Surface area as calculated for a delineated watershed basin for the area (Fig. 4) using the aforesaid discussed algorithms were used for pixel based comparative study.
The test area (for which a DEM was readily available) was chosen for its rough topography depicted by a greater range for slope (0 0 -79 0 ).
Results
Summarising the surface area value for the watershed basin we find the SHR method reflecting a higher value as compared to the slope method (Table 1).
The maximum area for a pixel using the slope method was much higher compared to the SHR method. This indicated that the SHR method was not sensitive to extreme slope values, which is natural owing to its dependency on neighbouring pixel values. Pixels with 0 slope value should ideally yield an SHR value of 1, which is actually not evident from the first record of Table 2. The reason is again attributed to the fact that SHR for 0 slope pixels were dependent on slope of neighbouring cells.
From our studies we could comprehensively conclude that SHR could be used to calculate surface area more accurately for a region with medium to low gradient. This finding makes it amenable that SHR method is better suited to regional and global scale studies whereas slope method can be applied to detailed scale (local) studies. SHR is important for terrain analysis where slope dependency between adjacent pixels is needed to preserve.
Discussion
GIS processing and decision making can be improved utilizing this surface area calculating functionality, especially for calculation of actual length and surface area for objects. Such functionality does add a new dimension to the accuracy of the surface area value although the dependency of accuracy also depends on the elevation model (TIN, DEM, DTED) in use, its information source (interpolated contours, stereo-imagery) and the resolution. Here we provide a comparison between two algorithms using terrain conditions that are amenable for both use. This comparison can also be extended to the study of intrinsic terrain roughness particularly for discriminating terrain types with multi-parameter terrain roughness estimates.
Slope can also be used to calculate the actual length of road and rail network as well, although a special algorithm can be devised to determine the departure of the object from one pixel to another. This departure (similar to aspect) can be used to calculate the apparent slope along that departure and from this apparent slope the actual length of the road or rail network can be approximated. This would enhance the present performance of area and length calculations in GIS and thereby enable decision making in applications involving cost associated area-length calculations.
References
- Burrough P.A. & McDonnell, R.A., 1998. Principles of geographical information systems: Spatial Information Systems and Geostatistics. OUP.
- Berry, J.K. 1987. Fundamental Operations in Computer-Assisted Map Analysis. Int. Journal of Geographical Information Systems, vol. 1, no. 2, pp. 119-136.
- Catmull, E.E., 1974, A Subdivision algorithm for Computer display of Curved Surface. PhD Thesis, University of Utah.
- Elghazali, M. S. & Hassan, M. M., 1986. A simplified terrain relief classification from DEM data using finite differences, Geo-Processing, v3 n2, 167-178.
- Horn, B.K.P., 1881, Hill shading and the Reflectance map, Proc. IEEE, Vol-69, No. 1 14-47.
- Onorati, G., Poscolieri, M., Ventura, R., Chiarini, V. & Crucilla, U., 1992. The digital elevation model of Italy for geomorphology and structural geology, Catena, v19 n2, 147-178.
- Skidmore A.K., 1989, Comparison of Techniques for Calculating Gradient and slope from a Gridded Digital Elevation Model. Int. Jour. Of G.I.S., v-3, No.4, 323-334.
- Strahler, A. N., 1952. Hypsometric (area-altitude) analysis of erosional topography, Geological Society American Bulletin, 63: 1117-1142.
- Travis, M.R., Elsner, G.H., Iverson, W.D. & Johnson, C.G., 1975. VIEWIT computation of seen areas, slope, aspect for land use planning. US Dept. of Agricultural Forest Service Technical report PSW 11/1975, Pacific Southwest Forest and Range Experimental Station, Berkley, California.
- Zevenbergen, L.W. & Thorne, C.R., 1987, Quantitative analysis of land surface topography. Earth Surface processes and Landforms, V-12, 47-56.S