Digital Mapping Techniques '03 — Workshop Proceedings
U.S. Geological Survey Open-File Report 03–471

Digital Karst Density Layer and Compilation of Mapped Karst Features in Pennsylvania

By Stuart O. Reese, P.G. and William E. Kochanov, P.G.

Pennsylvania Geological Survey, 3240 Schoolhouse Road, Middletown, PA 17057
Telephone (717) 702-2017; fax (717) 702-2065; e-mail: streese@state.pa.us, wkochanov@state.pa.us

INTRODUCTION

In 1985, the Pennsylvania Geological Survey began a series of investigations to map karst features of carbonate rocks in Pennsylvania (fig. 1). Locations of surface depressions, sinkholes, surface mines, and caves were compiled from municipal questionnaires, field surveys, published literature, and unpublished data sources as well as through an extensive aerial photograph review. Results of these investigations were released in a series of county-based open-file reports (Kochanov, 1987–1995). This information has been available since 1998 in paper maps and through an online database, but no digital layer files had been developed for geographic information system (GIS) tools.

Figure 1. Carbonate rocks in Pennsylvania (modified after Pennsylvania Geological Survey, 2000).

Over the course of several months, a digital compilation of karst data points was completed. In total, the compilation included over 111,000 data points from 14 counties and 107 7.5-minute quadrangles (fig. 2). Digital data allowed for GIS mapping and for computer analysis of karst features. A colorized density surface was created from these merged files using ArcView 3.2 software. Mapping the density of karst points is useful for assessment of potential structural and environmental problems associated with karst geology. High-density areas of karst points where land subsidence may be a problem are noted, or where karst features can serve as direct recharge zones to the groundwater. These areas are highly vulnerable to groundwater contamination.

Figure 2. Location of quadrangles with mapped karst features in south-central and southeastern Pennsylvania.

DIGITAL PROCEDURES

Data Compilation

Previously, karst feature locations from the open-file reports were plotted on a mylar stable base over the corresponding topographic map. The data points were digitized from the mylars using GSMAP v.8 (Selner and Taylor, 1992) and their coordinates were entered into a relational database management system. Karst data points were recorded in latitude/longitude(decimal degrees) by quadrangle. ArcView 3.2 was used to compile the karst point data from the database. ArcView shapefiles were created and merged to form a regional dataset. The shapefiles were placed over base maps of digital raster graphic images of 7.5-minute quadrangles in Universe Traverse Mercator (UTM), North American 1927 datum, Clarke 1866 spheroid, in UTM Zones 17 or 18.

The digital database was cross-checked against the original locations. As the data were reviewed, it became apparent that for many of the quadrangles, a systematic digitizing error had been introduced into the latitude and longitude data. These data points were corrected using the ArcView extension ShapeWarp (Version 2.2). In addition, on-screen digitizing procedures were used to create new files where previous digitizing had not been done. Compiled karst data points were identified by feature type (surface depression, sinkhole, surface mine, or cave), quadrangle, and county. County name was assigned by a spatial join command (ArcView geoprocessing, assign data by location).

Density Surface Preparation

The ArcView extension Spatial Analyst (version 2.0) was used to develop the digital density layer of the combined data points. A density surface is based on the division of the study area into square cells, which can be sized as appropriate. ArcView software calculates a density value for each cell by counting the number of points within a defined search radius from the center of each cell (fig. 3) and dividing by the search area. The density value (features per square mile) is assigned to the cell. The search circle is then shifted to the next cell and the floating process is repeated until all of the cells have been assigned a density value. This process smoothes the density layer over the study area.

Figure 3. Procedures used in density calculation.

ArcView has two options for density calculations — a simple density formula (described above) and a weighted “kernel” procedure. A weighted method called a kernel (a “quartic approximation to a Gaussian kernel”) can be used which assigns more value to points located closer to the center of the cell. For this project, the simple density function was used. The kernel method further smoothes the data but in this case caused a more bi-modal appearance of the density surface. For this reason, the simple density surface was retained.

Density calculations must be done in a map projection that minimizes error for area, distance, and direction. For example, a density calculation using a decimal degrees (geographic) framework would result in severe errors. The projection used for the density calculations was the Albers equal area conic projection with standard parallels at 40°N and 42°N, and a central meridian of -78°W. This projection maintains true area and shape with negligible distortion at the scale used, which was less than or equal to 1:24,000.

When doing density calculations, the recommended number of cells is between 10 and 100 cells per density unit (Mitchell, 1999). The density unit of features per square mile was used here, which, using this recommendation, would equal approximately 100 cells (with dimensions of 160 × 160 meters) per square mile. This results in a 5 MB file size with over 1.3 million cells. At a scale of 1:24,000 (the scale at which the features were mapped), the density surface has a noticeably blocky appearance. A smaller cell size will smooth the surface, but require more computer processing time and file storage space. A cell size of 25 meters was chosen to smooth the data to the lower limit of mappability at the 1:24,000 scale. This produced a 213 MB file with over 55 million cells. Despite the larger file size, the resulting density surface portrays a smooth gradation of cells and allows more local variation to be seen at the 1:24,000 scale.

The chosen search radius influences the appearance of the density surface. The larger the radius, the more generalized the patterns will be. The smaller the radius, the more local variation is portrayed, to the extreme of remapping the point data. The 250-meter search radius was selected in a trial and error process to show enough local variation without over-generalizing the density. Karst feature data points were overlaid upon the density surface to evaluate different parameter values. The 250-meter radius parameter provides for a smooth simple calculated density surface.

Further preparation of the density surface was accomplished using ArcGis 8.3. The cells were color coded using a graduated color scheme for the density values. The density color scale that was developed represents an ESRI ArcMap 8.3 “quantile gradation” of the density values. Quantile values are useful in comparing the density over an area, and from map to map. The large number of classes (30 quantiles) allows the color gradation to be displayed as a range and it allows the value of “0” to be given a transparent definition on the map. The mapped areas of higher density of karst features are portrayed in orange and red colors. Values in the red approach 640 karst features per square mile (one feature per acre) or more. The lowest density value, represented by thedarkest green color, indicates at least one karst feature within the 250-meter search radius of the cell (approximately 48 acres or 0.07 mi2).

The units of measurement were selected to consider the proposed audience of the product. Although units of meters were used to develop the density grid, density was calculated for each grid cell in units of karst features per square mile (1 mi² = 2.59 km²). A value in acres equivalent to the square mile value was added to the karst map. Units of square miles (and acres) allow the non-scientist to more easily relate to the map data. In the maps generated, yellow, orange and red colors approach 640 features per mile, or about one per acre.

DIGITAL PRODUCTS

The main products include digital coverages of karst points by county and by quadrangle, and a regional map showing the density of the karst features. Figure 4 shows an example of the density layer. County maps in a mapbook style were completed. Such easy-to-interpret maps will help homeowners, municipal planners and others understand the intensity of mapped karst features per square mile or acre. Because of potential misinterpretation, efforts were made to qualify the digital products.

Figure 4. Density of mapped karst features in the Lititz 7.5-minute quadrangle, Lancaster County, Pa.

CAVEATS AND CONCLUSIONS

Beyond the accompaniment of metadata with a GIS coverage, explicit caveat statements are very important because of the multiple meanings that could be interpreted from the karst data. On the regional karst density map, the bright colors are quickly noticed as areas of concern. However, the proper response needs to be guided. Banning all activity in the high-density zones is not a rational response, but neither should caution be thrown to the wind. The recommended uses of the data should be made obvious. Here, the karst data are useful for regional planning and preliminary site studies, but they are not a substitute for site-specific subsurface investigations.

The occurrence of a sinkhole, a subsidence feature that breaks the land surface, depends on numerous factors including rock type, geologic structure (the presence of fractures, joints and faults), soil cover, surface hydrology, and land use. Areas of subsidence are not necessarily restricted to the high-density areas shown in dark shades on the map. Surface depressions, which by definition do not show a land-surface break, were the dominant type of mapped karst feature (96 percent). However, subsidence can occur in areas where there are no discernible surface depressions, or where sinkholes are not observed.

In addition, there may be instances where subsidence features are shown outside the mapped limits of carbonate bedrock. Surficial material such as colluvium typically conceals the actual contact between non-carbonate and carbonate bedrock formations. Undetected faults also can displace carbonate bedrock and account for subsidence features outside the current mapped formational contacts of carbonate rocks.

Land use can bias the detection of the karst feature. Urban land cover often masks karst features, making them difficult to detect. Sinkholes, though highly visible and often disruptive when they occur, are typically quickly filled or covered. Therefore, mapped karst features are most often under-represented in the urban setting. Urban land cover often accelerates the formation of sinkholes through changes to surface water drainage. Land thought to be free of sinkholes may suddenly develop karst subsidence features, especially where the surface hydrology has been altered. In more rural areas, karst features can be difficult to discern in wooded areas, whereas karst features in fields are more easily detected.

All of the potential biases must be considered when using the karst density maps and the digital data of the karst features. Caveat statements are needed to provide proper direction on the use and interpretation of the digital products. Because the mapped features are based on interpretation, cautionary statements are extremely important to direct the use of such digital products, especially when there are known limits to the mapping process. An understanding of the limits of the data is crucial to the responsible use of the data.

ACKNOWLEDGMENTS

We wish to thank staff members of the Pennsylvania Geological Survey Dr. Jon Inners and Michael Moore for critically reviewing the manuscript.

REFERENCES

Kochanov, W. E., 1987a, Sinkholes and karst-related features of Lehigh County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 87–01, scale 1:24,000, 19 p, 6maps.

Kochanov, W. E., 1987b, Sinkholes and karst-related features of Northampton County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 87–02, scale 1:24,000, 24 p, 10 maps.

Kochanov, W. E., 1988a, Sinkholes and karst-related features of Berks County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 88–01, scale 1:24,000, 13 p., 16 maps.

Kochanov, W. E., 1988b, Sinkholes and karst-related features of Lebanon County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 88–02, scale 1:24,000, 11 p., 4 maps.

Kochanov, W. E., 1989a, Sinkholes and karst-related features of Cumberland County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 89–02, scale 1:24,000, 16 p., 9 maps.

Kochanov, W. E., 1989b, Sinkholes and karst-related features of Dauphin County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 89–01, scale 1:24,000, 8 p., 5 maps.

Kochanov, W. E., 1989c, Sinkholes and karst-related features of Franklin County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 89–03, scale 1:24,000, 15p., 20 maps.

Kochanov, W. E., 1990, Sinkholes and karst-related features of Lancaster County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 90–01, scale 1:24,000, 12 p., 18 maps.

Kochanov, W. E., 1993a, Sinkholes and karst-related features of Bucks County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 93–03, scale 1:24,000, 9 p., 4 maps.

Kochanov, W. E., 1993b, Sinkholes and karst-related features of Montgomery County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 93–02, scale 1:24,000, 7 p., 5 maps.

Kochanov, W. E., Lichtinger, J. F., and Becker, Mona, 1993, Sinkholes and karst-related features of Chester County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 93–01, scale 1:24,000, 9 p., 10 maps.

Kochanov, W. E., 1995, Sinkholes and karst-related features of York County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 95–06, scale 1:24,000, 9 p., 10 maps.

Kochanov, W. E., and Miller, Rachel, 1995, Sinkholes and karst-related features of Adams County, Pennsylvania: Pennsylvania Geological Survey, 4th ser., Open-File Report 95–05, scale 1:24,000, 9 p., 6 maps.

Mitchell, A., 1999, The ESRI Guide to GIS Analysis, Volume 1: Geographic Patterns & Relationships: Redlands, Ca., Environmental Systems Research Institute, Inc.

Pennsylvania Geological Survey, 2000, Limestone and dolomite distribution in Pennsylvania: Pennsylvania Geological Survey Map 15, scale 1:2,000,000. http://www.dcnr.state.pa.us/topogeo/maps/map15.pdf.

Pennsylvania Geological Survey, 2003, Density of Mapped Karst Features in South-Central and Southeastern Pennsylvania, by William E. Kochanov and Stuart O. Reese, 2003: Pennsylvania Geological Survey Map 68, scale 1:300,000. http://www.dcnr.state.pa.us/topogeo/map68.

Selner, G.I. and Taylor, R.B., 1992, GSMAP, GSMEDIT, GSMUTIL, GSPOST, GSDIG and other programs, version 8 for the IBM PC and compatible microcomputers to assist workers in the earth sciences: U.S. Geological Survey Open File Report 92–217.