Outcrop Prediction Problem 2 for QGIS 3.6



Outcrop Prediction Problem 2 for QGIS 3.6IntroductionOutcrop prediction is based on the concept of intersecting a planar surface (geological contact) with an irregular surface (topographic surface). In reality this occurs in 3D however the final depiction is on a 2D topographic map. The underlying assumption is that the geological contact or surface is a perfect plane. In many undeformed geological terranes this assumption is a valid over relatively small areas such as a 1:24K quadrangle. The value of the method is that it predicts on the topographic surface where the geological plane outcrops. This has obvious value in the case of a planar fault surface, or planar water table surface. If the method is applied to both the upper and lower planar contacts of a bed the result will be the predicted outcrop belt of the bed. The method also yields a structure contour map of the planar surface. This combined with the topographic map can be used to predict the depth below the ground surface of the plane at any point on a map – an obvious advantage when attempting to estimate costs of drilling wells or overburden removal in mining.Geometrical RelationshipsWhether the outcrop prediction is done “analog” with traditional orthographic projections, or “digital” with computer CAD/GIS, the fundamental process is the same – contouring a geological plane and then comparing those contours to the topographic contours. The structure contours of the geological plane will be straight lines, with the spacing between lines a function of the dip angle (inclination) of the plane:s = ci / tan(α)where s is the perpendicular distance between adjacent structure contours, α is the dip angle of the plane, and ci is the topographic contour interval. The elevation value of the structure contour always decreases by the contour interval value (ci) in the dip direction.If the outcrop of the structural plane occurs at a topographic elevation equal to a contour interval, the “s” spacing can be used from that point to construct all necessary structure contours as a set of parallel lines. If a contact can be traced on a map surface from aerial photography or from a geological map the “starting” structure contour can be set to where a topographic contour crosses the contact line. However, this may not be possible so the relationship:Tan (α) = (Δy)/(Δx) Can be used to calculate the exact offset from the control point to start the structure contour. For example, assume that the map topographic contour interval is 20 feet, and a bedding plane contact with attitude 330, 55 SW is discovered at 707 feet elevation. A 720 topographic contour is near the outcrop so the offset from the outcrop to the 720 structure contour on the contact is needed. First you should realize that the offset direction from the outcrop is in the 030 because that is perpendicular to strike and in the “up-dip” direction (i.e. elevation is gained from 707 to 720). The amount of map distance between the outcrop and 720 structure contour is calculated from:Tan(55) = 13/x x = 13/Tan(55) = 9.1 feetWhere x is the offset distance from the outcrop to the 720 structure contour. From this offset point the 720 structure contour can be drawn as a straight line striking 330. Where this line intersects the 720 topographic contour would generate outcrop control points. From the structure contour map the elevation of any point on the geological surface can be estimated, and this generates three possible scenarios when compared to the topographic contours/surface:The point is above the topographic surface (i.e. its elevation is higher than the topographic surface at that point) indicating that the geological plane has been eroded.The point is below the topographic surface (i.e. its elevation is lower than the topographic surface at that point) indication that the erosional surface is above the geological plane.The point on the geological plane is equal to the elevation of the topographic surface indicating the outcrop position of the plane. The line connecting these points is the trace of the outcrop on the map.Therefore, whether the process is accomplished manually or digitally finding the points where the structure contours intersect the matching topographic contours will trace out the outcrop pf the structural plane. If the outcrop position of the top and bottom of a tabular unit (i.e. top and bottom contacts are planar and parallel) are known, then the two sets of contours would allow tracing the contact of both planes across the topography. In this case the area between the two traces is the outcrop belt of the unit. Alternatively, if the outcrop position and orientation of the bottom of the unit where known, and the bedding thickness was measured and can be assumed to be constant, the structure contours of the top of the bed can be predicted by:w = t / sin(α)where w is the perpendicular distance between structure contours of equal value between the top and bottom contacts, t is the thickness of the bed, and α is the dip angle. The unit structural top contour would be offset “w” distance from the bottom contact equivalent contour in the dip direction. The spacing between adjacent top contact contours would be the same as calculated for the bottom contact because they both have the same dip and dip direction. Having both sets of structure contours overlaying the topographic contours can cause a map to become overly “complex” so whether the process is manual or digital it is advisable to use different color for each set. Figure 1 displays a topographic contour map of the USA campus at a 5-foot contour interval (black contour lines), and a structure contour map of the bottom of a planar clay unit (blue contour lines) using the same 5-foot contour interval. The dip angle of 1.43 degrees yields a spacing between adjacent structure contours of 200 feet. As you can see in Figure 1, the map can be “complex” even with only one set of structure contours superimposed on the topographic contours. As you can imagine to process of tracking where structure contours intersect the same elevation topographic contour is tedious, time consuming, and error-prone. In other words, a perfect task for a computer!Using QGIS to Process Outcrop PredictionThe good news about outcrop prediction is this- QGIS can easily solve the outcrop prediction. The steps for calculating the outcrop of a planar geological surface with ArcGIS are summarized below:Step 1: Add base map feature to QGIS project.Step 2: Use survey data to interpolate an elevation grid that covers the area of interest, or you may be able to download a DEM from online sources.Step 3: Use the Python programs “PlaneTrendSurf.py” to generate an elevation grid of the geological plane based on its outcrop position and orientation. If you need to solve a 3-point problem do this now because you will need the trend and plunge of the true dip as input into this program. Do the same for the other contact if modelling a tabular unit.Step 4: Use the “Raster Calculator” to subtract the topographic grid from the planar surface grid to produce “residual” rasters for top and bottom contacts.Step 5: Reclassify both of the of the residual raster features produced in step 4 into simplified raster features where the pixel value “1” represents strata above the contact plane (i.e. negative residuals), and “2” represents strata below the contact (i.e. positive values). (Processing Tools: SAGA > Raster Tools > Reclassify values).Step 6: Use “Raster Calculator” to add together the two rasters produced in step 5. The composite will contain three possible pixel values: 2=strata above top contact, 3=strata between bottom and top contact, and 4=strata below bottom contact.Step 7: Use the SAGA routines “vector <-> raster > contour lines” and “vector line tools > smooth lines” to generate lithologic “contacts from the geology raster created in step (6). Use the methods covered in the example problem to create the lithologic polygons with “area”, “perimeter”, and “lithocode” attributes.Step 8: Use the attribute table “field calculator” in the “Lithology” polygon layer to fill in the values of “area”, “perimeter”, and “lithocode” as in the example problem.Problem 2: Outcrop Prediction of a Clay Unit on USA CampusIn this problem you will be given a campus base map containing roads, buildings, etc., along with topographic survey points and two locations where the bottom and top of a tabular clay unit is exposed. The strike and dip measurements on each contact are embedded in the attribute table of the given points. Drilling on campus has verified that a sandy unit underlies the clay, and a gravel unit is on top of the clay. The first task is to build a geologic map of the campus area using QGIS outcrop prediction methods for the three units, and then answer the questions at the end of this discussion using QGIS tools. At the start of the project you will be given a digital elevation model (DEM) raster named “Elevation” that was derived from the supplied topographic survey.Step 1: Add Base Map Features to New ArcMap Project.Download the following files from the GIT 461/561 USA online site from the “Resources” section:OP_prob2_start.zip {starting files for problem 2}PlaneTrendSurf.py{Python script for calculating a planar trend surface grid}Alternatively, if you don’t have a USA online account you can access all of the starting files at the below web site: for the “Outcrop Prediction Problem 2 QGIS 2.18” Heading in the list of projects.The base map features below will be made available to you to start the project in the form of the below .shp or .tif files that are unzipped from the archive. Make sure that you create a folder and unzip all of the files from the ZIP archive to the folder:OP_prob2_start.zip >>Buildings.shp (campus buildings)Mine.shp (proposed gravel mine area)Roads.shp (campus road system)Stations.shp (outcrop points for bottom and top of clay)SurveyPts.shp (topographic survey points used to build DEM)Topocontours.shp (topographic contours generated from survey points)Elevation.tif (elevation DEM generated from survey points)Border.shp (mapping boundary)You will also need access to the below Python program:PlaneTrendSurf.py {download directly from the web site}Start QGIS and create a new project file that uses a CRS of AL West state plane coordinates (feet) in NAD27 map datum. Add the base map features listed above to the project file. Figure 2 displays the QGIS project file with these features added to the project. The “Identify feature” tool window is open to display the data attached to the NW outcrop point. The outcrop in the NW corner of the map is the exposed bottom contact of the clay unit. The outcrop in the NE portion of the map is the exposed upper clay contact.The elevation raster was generated with the “SAGA > Raster Creation > Multilevel B-spline interpolation”. The elevation grid was constructed from the “Z” value field in “Survey Pt.s” using the following settings:Grid spacing: 25Extents: SW corner = 287414.5, 253440.5 ; NE corner = 287714.5, 253815.5The “Topocontours” vector layer was generated from this elevation grid with the SAGA “contour lines” procedure. Step 2: Generate a Plane Elevation Grid with PlaneTrendSurf.Start the “Python” console from the QGIS toolbar and load the “PlaneTrendSurf.py” program into the script editor with “Load Script”. Scroll down through the code to check the “path”, “ofn”, and “rfn” variables to make sure they point to the proper input (rfn) and output (ofn) file names. An example input file is loaded into NotePad in Figure 3 along with the “PlaneTrendSurf.py” program displayed in the QGIS Python console window. Note the data file format:# Line 1. trend and plunge degrees of the true dip vector (trend in azimuth format). (ex. 45.0, 66.4).{find attitude from stations layer attribute data}# Line 2. Elevation (Z value) of the measured attitude. (ex. 1200.0).For consistency use:190 (bottom)150 (top)# Line 3. "X,Y" coordinates of the position of the attitude measurement (ex. "300000.0 , 4000100.0").For consistency use:283552.5, 256783.8 (bottom) 288233.0, 255241.5 (top)# Line 4. SW corner (lower left) "X,Y" coordinates of the calculated raster image. (ex. 280000.0, 4000000.0).NOTE: for consistency with the key use these coordinates for the SW corner point: 282889.5, 251540.5# Line 5. "Columns,Rows" in calulated raster (ex. "100, 100").Note: for consistency with key use :217, 218# Line 6. Grid spacing between rows and columns of raster grid (ex. "50.0").For consistency use:25.0Use the “Identify” tool to determine the coordinates of the two given data points. The attitude field of the “Outcrop stations” feature class contains the attitude in trend and plunge of true dip vector. The above example values DO NOT MATCH your data- determine the proper values with the “Identify” tool or by opening the attribute table for “Stations”. You can create the input text file with any text editor, such as “Notepad”.When the input text file is ready (in this example “prob2_bot_in.txt”) use the “run” menu in Python to execute the “PlaneTrendSurf.py” program. The last screen of results from the program are the coefficients of the 1st order (planar) trend surface: Z = c0 + c1*x + c2*yWhere z is the elevation, and x and y are the coordinates of the map. The output grid text file (“prob2_bot_grd.txt”) is the grid produced by using this equation to calculate the elevation at each grid node.The next portion of this step is to import the output grid – simply use the “add raster layer” tool to add this grid as a new layer. Figure 4 displays the appearance of the clay bottom or top contact grid- they both look the same because they are parallel planes. However, you should check to make sure they are consistent with the original data control points. Note that the strike of the unit is northeast parallel to bands of color, with color values representing decreasing elevation (darker) to the southeast dip direction. Step 3: Subtract the Topographic grid from the Clay Bottom and Top Contacts to create Residuals.Use “Raster > Raster Calculator” from the main menu to subtract the elevation grid from the clay bottom planar grid. The results should be added to your project as “resid_bot” and “resid_top” respectively. Check each grid with the “Identify features” tool to ensure that near the original control points the residuals are near 0.0 in value.Step 4: Reclassify the Two Residual Rasters.To reclassify the two residual grids find the “processing tools: SAGA > Raster tools > Reclassify values” tool and start it. Reclassify the bottom clay residual grid as specified in Figure 5. Use the same procedure to reclassify the top clay outcrop grid replacing “1” for negative values (strata below contact outcrops) and “2” for positive values (strata above contact outcrops). The new reclassified raster will be automatically added to the project. Note that these two rasters display a lithologic outcrop map – values of “1” above the contact plane and values of “2” below it. The “Reclass_Top” should appear as in Figure 6.Step 5: Add together the Top and Bottom Reclassified GridsUse “Raster > Raster Calculator” to add the top and bottom contact reclassified grids together to produce a composite grid (“GeologyRaster”). In the composite grid “GeologyRaster” pixel values of “2” will represent the outcropping of the upper gravel, the value “3” the exposure of the clay unit, and “4” the lower sand unit. Your project map should appear similar to Figure 7. Step 6 Convert the Raster Composite to a true Polygon Topology Use the same procedure that was outlined in problem 1 and the example problem to convert the “GeologyRaster” to a true polygon layer:Use “processing tools: SAGA > vector <-> raster > contour lines” to generate contact contours.Use “processing tools: SAGA > vector tools > smooth lines” to smooth the contours. Name this layer “contacts”.Use editing and snap modes to complete the contact line topology. Make a new “Lithology” polygon layer that covers the mapping area. Use “contacts” to cut the “Lithology” layer into geologic polygons.Use “Field Calculator” to fill in the area and perimeter values. Use the attribute table in edit mode to classify the polygons as “sand”, “clay” and “gravel”.The polygon topology will enable the ability to process calculations involving the area of the polygons. Open the attribute table for “Lithology” and note that “area” has already been added as a field. In edit mode add another field named “Lithology” as a text field of length 16. Use the layer symbology to color the “Lithology” layer:Upper Gravel strata: greenMiddle clay strata: redLower sand strata: yellowYour final lithologic polygons should appear similar to Figure 8.Step 7: Setup the Layout for Printing As in problem 1 use the layout composer to generate a layout with a title, legend, scale bar, and north arrow at a scale of 1:10,000. Print a copy to turn in with the below problem answers. Figure 9 displays the final layout composition.Questions Section2.1: What is the area of the clay unit outcrop belt in m2 ? __________________2.2: What is the thickness of the clay unit in feet? ________________2.3: To what depth would a driller have to drill to encounter the top of the clay unit if a vertical drill hole was started at the NW corner of the Administration Building? ___________________2.4: USA wants to mine the gravel unit for future construction materials. Add the “Mine” feature to your project to display the area to be mined. Calculate the volume of material in cubic feet:Cubic feet of gravel = __________________________________________{Hint: Use an approach to the one in the example problem to calculate this quantity. You can “clip” both the elevation and top contact grids to within the mine area. Subtracting the mine top contact from the mine elevation surface would yield a “Mine_Diff” raster that has pixel values equal to the difference (depth) to the top contact at the x,y position of the pixel. If you then multiply: (area of the mine) * (average pixel value of “Mine_Diff”) = volume of gravel To clip a raster layer with a polygon layer use “Processing Toolbox: SAGA > Vector - Raster > Clip Raster with Polygon”. 152400128270Figure 1: Problem 2 USA campus topographic map with manual structure contours superimposed.Figure 1: Problem 2 USA campus topographic map with manual structure contours superimposed.10728962548128005840730Figure 2: Starting files loaded into project and NW outcrop point selected with "Identify features" tool.Figure 2: Starting files loaded into project and NW outcrop point selected with "Identify features" tool.centercenter005840730Figure 3: Python program loaded into QGIS console with clay bottom input file displayed in Notepad.Figure 3: Python program loaded into QGIS console with clay bottom input file displayed in Notepad.centercenter005780405Figure 5: Appearance of top or bottom contact elevation grid.Figure SEQ Figure \* ARABIC 5: Appearance of top or bottom contact elevation grid.centercenter03048007308215Figure SEQ Figure \* ARABIC 6: Reclassify dialog window for the top contact residuals.Figure 6: Reclassify dialog window for the top contact residuals.centercenter0005133975Figure 7: Appearance of the "Reclass_top" raster.Figure SEQ Figure \* ARABIC 7: Appearance of the "Reclass_top" raster.centercenter005780405Figure SEQ Figure \* ARABIC 8: Composite raster "GeologyRaster".Figure 8: Composite raster "GeologyRaster".centercenter005780405Figure SEQ Figure \* ARABIC 9: Final symbology on "Lithology" polygon layer.Figure 9: Final symbology on "Lithology" polygon layer.centercenter006237605Figure SEQ Figure \* ARABIC 10: Layout composition for problem 2.Figure 10: Layout composition for problem 2.centercenter0 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download