Our project to detect relict charcoal hearths (RCHs) relies on using LiDAR data that shows the what remains of these hearths on the surface of the land. These hearths were constructed in forested areas near the trees that were felled to build them. As forests regrew, remains of RCHs were often preserved. The distinctive flat, 10-15m circular areas are visible in slope images derived from LiDAR data. Not all of these objects in forested areas are RCHs, but in Pennsylvania, many are. The density of RCHs in parts of the state is due to the extent of historical charcoal making to fuel the early steel industry. This was before steel making transitioned to using mineral coal starting in the mid-nineteenth century and completing in 1945.
In non-forested areas it is much less likely that similar looking objects are RCHs. The surfaces of non-forested areas today are more likely to be disturbed by development or agriculture since the era of charcoal making. Thus, any evidence of an RCH, even if once present, is likely to have been plowed or bulldozed over.
Given this premise, when detecting RCH it is more efficient to spend effort only on searching forested areas. At the same time, by ignoring detected objects that look like RCHs in non-forested areas, high likelihood false positive RCH detections can be removed from consideration. For these reasons, I wanted to have a vector layer of forested areas of Pennsylvania. Creating the vector layer takes several steps, described below.
PAMAP Program Land Cover for Pennsylvania, 2005
I based the forest vector layer on the PAMAP Program Land Cover for Pennsylvania, 2005. https://www.pasda.psu.edu/uci/DataSummary.aspx?dataset=1100. This is a detailed, large file size raster that describes the land cover of the state. Due to the detail, converting this entire raster file into polygons of a vector file consumed too many resources for QGIS to complete this on my machine. To solve that problem, I clipped the dataset into a smaller geographic area and then filtered the results to just forest land cover.
Clip by extent
Our project is conducting a search of RCHs in the entire state of Pennsylvania. To break the project down into smaller steps, the project searches the state by Forest District. Presently, we are searching Forest District 17 in the south east part of the state. Below is a screen shot of the land cover raster in gray scale and Forest District 17 outlined in green. The Forest District file is filtered on “DistrictNu” = 17 in QGIS.
Below is the detail of the QGIS command used to clip the raster. In QGIS’s menu it is: Raster | Extraction | Clip Raster by Mask Layer…
Input layer: ‘E:/a_new_orgs/carleton/pennsylvania_michaux/land_cover/palulc_05_utm18_nad83/palulc_05.tif’
Mask layer: ‘E:/a_new_orgs/carleton/pennsylvania_michaux/DCNR_BOF_Bndry_SFM201703/DCNR_BOF_Bndry_SFM201703.shp|subset=\”DistrictNu\” = 17’
Algorithm 'Clip raster by mask layer' starting…
Input parameters:
{ 'ALPHA_BAND' : False, 'CROP_TO_CUTLINE' : True, 'DATA_TYPE' : 0, 'EXTRA' : '', 'INPUT' : 'E:/a_new_orgs/carleton/pennsylvania_michaux/land_cover/palulc_05_utm18_nad83/palulc_05.tif', 'KEEP_RESOLUTION' : False, 'MASK' : 'E:/a_new_orgs/carleton/pennsylvania_michaux/DCNR_BOF_Bndry_SFM201703/DCNR_BOF_Bndry_SFM201703.shp|subset=\"DistrictNu\" = 17', 'MULTITHREADING' : False, 'NODATA' : None, 'OPTIONS' : '', 'OUTPUT' : 'TEMPORARY_OUTPUT', 'SET_RESOLUTION' : False, 'SOURCE_CRS' : QgsCoordinateReferenceSystem('EPSG:26918'), 'TARGET_CRS' : None, 'X_RESOLUTION' : None, 'Y_RESOLUTION' : None }
GDAL command:
gdalwarp -s_srs EPSG:26918 -of GTiff -cutline C:/Users/User/AppData/Local/Temp/processing_gEEUVX/081df2395865469b973df34e0f45f1e5/MASK.shp -cl MASK -crop_to_cutline E:\a_new_orgs\carleton\pennsylvania_michaux\land_cover\palulc_05_utm18_nad83\palulc_05.tif C:/Users/User/AppData/Local/Temp/processing_gEEUVX/899b60ac78354e019d9c2a2992ae4905/OUTPUT.tif
GDAL command output:
Copying raster attribute table from E:\a_new_orgs\carleton\pennsylvania_michaux\land_cover\palulc_05_utm18_nad83\palulc_05.tif to new file.
Creating output file that is 5681P x 4602L.
Processing E:\a_new_orgs\carleton\pennsylvania_michaux\land_cover\palulc_05_utm18_nad83\palulc_05.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.
Execution completed in 4.94 seconds
Results:
{'OUTPUT': 'C:/Users/User/AppData/Local/Temp/processing_gEEUVX/899b60ac78354e019d9c2a2992ae4905/OUTPUT.tif'}
Loading resulting layers
Algorithm 'Clip raster by mask layer' finished
To filter just the forest values, I referenced https://www.pasda.psu.edu/uci/FullMetadataDisplay.aspx?file=palanduse05utm18nad83.xml which shows values for
41 – Deciduous Forest
42 – Evergreen Forest
43 – Mixed Deciduous and Evergreen
QGIS’ raster calculator is used to the value of raster pixels to 1 if the land cover value is >= 41 and <= 43 (forested) using this formula:
1*("Clipped (mask)@1" >= 41 and "Clipped (mask)@1" <= 43)
Sieving out small areas of forest.
For the purposes of this project searching small areas of forest is unlikely to result in RCHs, thus searching these small areas may be a waste of time. To reduce that waste, the small forested areas are sieved from the raster.
Raster | Analysis |Sieve… Threshold of 9.
This raster has values 1 for forested and 0 for not. I only want to work with the forested data and so I wish to remove the 0 values. To do this in QGIS, right-click on the layer: Export | Save as… and specify 0,0 as No data values.
Create the vector layer
To create the vector layer from a raster in QGIS: Raster | Conversion | Polygonize… Name of field to create: “forested” gdal_polygonize.bat E:\a_new_orgs\carleton\pennsylvania_michaux\land_cover\d17_forest_cover_9july_just_1.tif C:/Users/User/AppData/Local/Temp/processing_exTjDd/2e62017e5fa3420b8dbfff6a5508356f/OUTPUT.gpkg -b 1 -f “GPKG” OUTPUT Forested
The conversion of a raster to vector layer takes a few minutes to process. However when I used the whole raster, before clipping, merging values, sieving and setting 0 = no data, making a vector layer took more than 1 hour and failed.
