QGIS Outcrop Prediction Example Problem



QGIS Outcrop Prediction Example ProblemIntroductionThis document discuss the step-by-step example of using QGIS with certain “plugin” extensions to process a outcrop prediction geologic map starting with a digital elevation raster and the orientation of the bottom and top contacts of a inclined tabular geologic unit. A Python application will create the top contact planar elevation raster file based on the attitude of the top contact and a single geographic position on the map where the attitude was observed. The bottom contact raster is derived from to top contact raster based on the thickness of the unit.Setup and Starting Data You can download the starting data from the below web site: for the section titled “Outcrop Prediction” on the web site and then download from the below links:OPexample_QGIS_start.zip {starting data files}PlaneTrendSurf.py {python program}Create a folder on a disk and save these files to that working folder. It is advisable to put any files created during this exercise in that same folder. If you are working on a workstation on campus note that most don’t allow the permanent storage of student files on the hard disk so it is advisable to create a folder on your flash drive and work in that folder periodically saving results to that folder.Getting Started Open QGIS form the desktop or start menu. Select “Project > New” from the main menu to create a blank map. Set project coordinate system using the “Project > Properties > CRS“ menu selection. Set the coordinate system to UTM NAD27 zone 16. Load the following files using either the “Layer > Add layer > add vector layer”:SurveyPts {points}Border {line}Bedding {points}Exprob_MineArea {polygon}After adding the above layers you should see the map border, a number of survey points, and one bedding data point. Right-click on the “SurveyPts” layer name, select “Open Attribute Table” and view the organization of the survey point data. The x and y coordinates are UTM zone 16 in meters, and the “Z_coord” field is elevation in meters. Close the attribute table and then open the attribute table for the “bedding” layer. Note that in addition of the x,y,z coordinates of the data point, the azimuth of the strike, the dip amount and dip direction for bedding at that point are entered in the table. This point represents the top of the stratigraphic unit to be mapped.Generating the Elevation GridAn elevation grid is needed to calculate the outcrop prediction so as a 1st step we will generate a grid from the randomly located survey data. To generate the grid select from the main menu “Processing > Toolbox” which should open a “Processing Toolbox” window (be aware that it may already be open). Find the “SAGA” modules and expand the “Raster Creation Tools” group. Find and run the “Multilevel b-spline interpolation” module and fill in the dialog window as in Figure 1. Note that the grid cell size is set to 1.0. After successfully running this routine an elevation grid will be created and added as a layer to your project. Rename this grid to “Elevation”. Generating the Top and Bottom Contact Elevation GridsThe python program “PlaneTrendSurf.py” will be used to generate the planar top and bottom contact grids. The input file for the program is described below:---------------------------------------------------------------------------------------------------------Line 1. trend and plunge degrees of the true dip vector (trend in azimuth format). (for this example use “90.0, 40.0”).Line 2. Elevation (Z value) of the measured attitude. (for this example use “110.0”).Line 3. "X,Y" coordinates of the position of the attitude measurement (for this example use "500165.0, 3000085.0").Line 4. SW corner (lower left) "X,Y" coordinates of the calculated raster image. (for this example use “500000.0, 3000000.0")Line 5. "Columns,Rows" in calulated raster. (for this example use "195, 145").Line 6. Grid spacing between rows and columns of raster grid. (for this example use "1.0").File name: ProbEx_top_in.txt. {below is input file layout:}90.0, 40.0110.0500165.0, 3000085.0500000.0, 3000000.0195, 1451.0{text file output after running “PlaneTrendSurf.py” with the above input file:}Output: ESRI raster file with following structure:Line 1: "ncols=" {number of columns in raster grid}Line 2: "nrows=" {number of rows in raster grid}Line 3: "xllcorner=" {x coordinate of the lower left (SW) corner of the grid}Line 4: "yllcorner=" {y coordinate of the lower left (SW) corner of the grid}Line 5: "cellsize=" {cell size spacing between columns and rows}Line 6: "nodata value=" {value that represents no data at grid node}Line 7: 1st z valueLine 8: 2nd z valueLine (ncols * nrows + 6) : (ncols * nrows) z valueOutput file name: “ProbEx_top_grd.txt”You can use “Notepad” or any similar text editor to create the input file. Just remember that comas are used to separate two values that occur on the same line, as in line 3 that contains the x and y coordinates of the outcrop point. This means that you cannot use comas within the x or y coordinate number.The input file for the bottom contact – “ProbEx_bot_in.txt” in this example – will be identical to the top contact input file except that the elevation will be a lower value based on the thickness of the tabular unit:d = t / cos(δ)where d represents the decrease in elevation, t is the thickness of the tabular unit, and δ is the dip angle. For this example the thickness is 12.75m therefore:d = 12.75 / cos(40) = 16.6mThe “ProbEx_bot_in.txt” file should appear as below:90.0, 40.093.4 {i.e. 110m – 16.6m = 93.4m}500165.0, 3000085.0500000.0, 3000000.0195, 1451.0Note that the only difference between the input files is the elevation of the control point on line 2. Proceed to create the ASCII elevation rasters “ProbEx_top_grd.txt” and “ProbEx_bot_grd.txt” from the above input files.From within QGIS use “Layer > Add Layer > Raster Layer” to add the “ProbEx_top_grd.txt” and “ProbEx_bot_grd.txt” files as 2 layers. Both of the grids will appear the same because they are parallel planar surfaces. The top contact surface grid is displayed in Figure 2. Note that the “Identify Results” windows shows a “click” near the bedding data point on the top contact surface grid, and it is near 110 elevation as it should be. Also note the pattern of the grid – decreasing in elevation from west to east along constant elevation color bands running north-south parallel to strike. At this point it would be advisable to check the bottom contact before proceeding to the next step. The bottom contact should be 16.6m below an equivalent (x,y) point on the top contact.Create the Residuals by Subtracting the Elevation from Top and Bottom GridsResiduals are calculated by subtracting the elevation grid from the top or bottom grid to produce “top_resid” and “bot_resid” raster grid feature classes respectively. Both of these rasters have the following characteristics:Grid NodeGeometry> 0.0Contact is above topographic surface therefore at this grid point the strata below the contact is exposed at the ground surface.= 0.0Contact outcrops at this point.< 0.0Contact is below topographic surface therefore the strata above the contact at this grid point is exposed on the ground surface.Note that if you symbolized all values < 0.0 with a blue color, and values > 0.0 with red, and 0.0 pixels black you would have a basic geologic map with blue representing the strata above the contact outcropping at the surface, red would be the outcrop of the strata below the contact, and the black pixels would trace the contact itself. At this time run the “Raster > Raster Calculator” selection that will display the window dialog in Figure 3. Fill in the dialog as indicated to produce the “top_resid” top contact – elevation residuals grid.Figure 4 displays the top contact grid residual “top_resid” in a gray-scale color scheme. Note that in the “Identify Results” window a left-click near the bedding data point yields a residual near zero at the outcrop point.Reclassify the Residuals The residual rasters will be reclassified so that values <= 0.0 are assigned a value of 1, and values > 0.0 are assigned a value of 2. After this step both “top_reclass” and “bot_reclass” will have the following characteristics: SEQ CHAPTER \h \r 1Grid Node ValueReclassified Value< 0.0 1 (strata above contact exposed)>= 0.02 (strata below contact exposed)Find the SAGA “Raster Tools > Reclassify” tool in the Processing menu or Processing Toolbox window. Figure 5 contains the window dialog settings for reclassifying the “top_resid” to make the “top_reclass”. Make sure that you set the method to “Range”. For minimum value you just need to set a value that is guaranteed to be below the minimum residual value. The “top_reclass” map is displayed in Figure 6. Note that it consists of only two values, 1 for the outcrop of strata above the top contact surface, and 2 for strata below the contact. Add the Two Reclass RastersThis simplification will allow the top and bottom rasters to be “added” together to produce a raster that has the below characteristics: SEQ CHAPTER \h \r 1Top Grid+Bottom Grid=Composite Grid112=younger strata exposed(strata above top and bottom contact)213=intermediate strata exposed(strata below top contact, above bottom contact)224=older strata exposed(strata below top and bottom contact)Note that the condition of the top grid node = 1 and the bottom grid = 2 at the same map position is not possible if the top contact is structurally above the bottom contact. The composite grid separates the older, intermediate, and upper stratigraphic units into three unique values, therefore, it represents the geologic map in raster form. Use the raster math dialog to add the “top_reclass” and “bot_reclass” to produce “Litho_raster”. Figure 7 displays this composite as the raster geologic map. Note that the lower value pixels are stratigraphically younger therefore the colors in the “layer s panel” window are in proper stratigraphic order.Converting the Raster Geologic Map to Lithologic PolygonAlthough the “Litho_raster” is a raster geologic map there are several reasons why you may want to convert the raster into a polygon feature class. For example, if converted to a polygon feature the area of each polygon entity is automatically calculated whereas a fairly complex analysis would have to be processed to yield the same value with a raster. The boundaries between polygons are smoother and give a better quality to plotted maps. The first step in producing a true polygon topology from the “Litho_raster” layer is to generate contact lines by “contouring” the raster values contained in “Litho_raster. Because this raster contains 2, 3, and 4 as raster values we will generate contours at 2.5 and 3.5. Start the SAGA command “Vector <-> Raster > contour”. The window dialog for this procedure is displayed in Figure 8 – fill in values as indicated. The contour values generated will be jagged because they follow the raster pixel boundaries so the contours need to be “generalized” – run the SAGA procedure “vector line tools > simplify line” using the “contours” as input. The window dialog is displayed in Figure 9 – the output smoothed lines are named “contacts”.The next step constructs a “lithology” polygon that initially covers the entire map area, but is “sliced” by the contact lines into 3 separate polygon parcels. In this example the polygons from left-to-right contain the lower (value=4), intermediate (value=3), and upper (value=2) stratigraphic units. The construction steps are in order:Create the initial polygon rectangle (SW corner= 500000,3000000; NE corner =500195,3000145).Use the shape edit toolbar function “cut polygon with lines” to slice the polygon with the “contact” layer lines.Add area, perimeter, and lithology type attributes to the “Lithology” polygon layer.From the main menu use “Layer > New Layer > New shape file layer” to open the layer creation dialog. In this dialog you should indicate that you are creating a new polygon layer with attributes “area”, “perimeter”, and “lithocode”. Area and perimeter are decimal numbers; “lithocode” is a text attribute. Fill in the dialog as indicated in Figure 10. At this point make sure the new “lithology” polygon layer is highlighted in the layer window, and then select the “editing” tool from the digitizing toolbar. As always, you may have to “turn on” specific toolbars by right-clicking on the unoccupied area of the toolbar. Selecting this tool will toggle “on” editing of features in this layer, and make it possible to create new items. Turn on the “Snapping” toolbar at this time and make sure that snapping is “on” and that snapping to “vertices” is active. This will make left-clicks to automatically “snap” to vertices of lines. Also make sure the “Border” layer is on because you will be snapping to the corner vertices. Click on the “Create Polygons” tool in the digitizing toolbar. Proceed to snap to the 4 corners and then right-click to finish. You can leave the attribute fields blank (or NULL) for now. You now need to slice the starting polygon using the contact lines as cutting edges. Select the following objects taking the below actions:Highlight the “contacts” layer in the layer list window.Use the “selection” tool to highlight one of the contact lines. The line should turn yellow (selection color).Highlight the “lithology” polygons layer in the layer list window.Use the “selection” tool to highlight the lithology polygon. It should turn yellow – this will make the selected contact line disappear but no worries.Select the “Cut features using a line from another layer” button in the “Shape Digitizing” toolbar. This should slice the polygon into two portions divided by the contact line that was selected. You will have to indicate the layer where the highlighted line exists – in this case “contacts”.Use the selection tool to highlight the west “slice”. It should look similar to Figure 11. Proceed to use the other contact line to create the 3 lithologic polygons. At this point it would be nice to assign different colors to each polygon (younger, middle, older strata) based on a text field value. Right-click on the “Lithology” name in the “Layers” window, and then select “Attribute Table”. The attribute table will now display – highlight the row in the attribute table and then fill in the “Lithocode” values for “Older” (west), “Intermediate” (middle), “younger” (east). Now right-click on the “Lithology” layer name again and select “properties”. Then select the “Style” tab. Figure 12 contains the window dialog filled out for this layer style – see if you can fill it out the same way. Work from top down being sure to select “Categories” and the “Lithocode” column. Then select “Classify” to list all unique values in the “Lithology” field. You can double-click on the color squares to change them to specific colors. Now find the field calculation query tool in the attribute table (activate edit mode if necessary) and start the calculation query dialog. Look for the geometry variable “$area” to add to the query and then click “OK” to run. The field will automatically fill in with the area values in square meters. Figure 13 contains the calculation query dialog window. Do the calculation for the “perimeter” attribute also. Be aware that if your attribute table has row(s) selected, the query calculator will only apply the calculation to the selected rows. You can quickly un-select any rows with the “deselect selections” button in the attribute posing the Final Map ProductQGIS can easily add the final touches to the geologic map. First, use the “Project > New Print Composer” to name a layout and then open a blank layout window. Use the “Layout > Add Map” to add a map frame to the layout. Proceed from this same menu to add:North Arrow {Add an arrow and then a large “N” label at the base of the arrow}Scale Bar{Set a scale of 1:1000}Title {Problem ? Geologic Map}Legend{Include only the vector layers- see Figure 14}You should wind up with something similar to Figure 14.CalculationsWith the lithologic polygon feature you can do a variety of simple calculations:Problem 1: What is the map area of the:Younger stratigraphic unit? 3792.4 m2 intermediate stratigraphic unit? 5732.7 m2older stratigraphic unit? 18998.0 m2Solution: Use the “Identify” tool to click on the unit polygon. When the results window is displayed look for the value of “area”. The area units will always be the square of the unit used by the map coordinate system – in this case UTM that uses meters. Recall that the “area” field of the Lithology polygon layer was added in a previous step, and was set to the value of “$area” with a calculation query.Problem 2: Assuming the intermediate unit is composed of coal what volume in m3 could be extracted from the map area? 13,200.3 m3The volume of material between the top and bottom contacts of the intermediate unit can be calculated by first estimating the map area of the mine, and then multiplying that value by the vertical depth between the top and bottom contact. In effect multiplying the map area by the thickness of the unit.Step 1: Find the area of the mine polygon by right-clicking on the “Exprob_MineArea” layer and selecting “attribute table”. The area of the single polygon will be one of the field values in the table (795.2m2). Step 2: Calculation of the depth from the top to the bottom contact is easily done from the equations for the planar trend surfaces output by “PlaneTrend.py”:Bottom: C0 = 419781.7; C1 = -0.839; C2 = -1.490Top: C0 = 419798.3; C1 = -0.839; C2 = -1.490Notable is the fact that C1 and C2 are the same for both planar contacts but this must be true if both contacts are parallel. The only difference in the 2 equations is the Z intercept that differs by a value of 16.6m. This is the vertical depth from the top to the bottom contact, therefore the volume is:795.2 m2 * 16.6 = 13,200.3 m3Problem 3: If the younger stratigraphic unit consisted of sand and gravel how much material (m3) could be mined from the mine area? 9702 m3Step 1: Use the “Raster > Extraction > Clip raster by mask layer” to clip out the section of the elevation grid that corresponds to the mine area. In the dialog use these parameters:Input file: ElevationOutput file: Mine_ElevationMask layer: Exprob_MineArea(check the “crop to extent of the crop line}This will extract the portion of the elevation grid that corresponds to the mine area on the map. Use the same extraction process to produce Mine_Top, the portion of the top contact that falls in the mine area.Step 2: Use the “Raster >Raster Calculator” to calculate the following:Diff1 = Mine_Elevation – Mine_TopMake sure that you use band1 (@1) in both cases. Also set the extents of the “Diff1” result to the “Current layer extent”. You should realize that the result of this calculation represents the depth from the surface down to the planar bottom of the younger stratigraphic unit. Step 3: Right-click on the “Diff1” layer name in the layer window and select “properties”. Scan the list to find the x and y dimensions of the raster. Also find the “pixel” area size. The map area of the mine equals:X dimension * y dimension * pixel area If the above value is multiplied by the average value of the pixel depth values the volume of the sand and gravel will result. Scan the metadata list for the average value of the raster. The final calculation is:18 * 44 * 1m * 12.25m = 9702 m307221220Figure 1: Multilevel B-spline Interpolation window dialog.Figure SEQ Figure \* ARABIC 1: Multilevel B-spline Interpolation window dialog.centercenter050805780405Figure 2: Appearance of the top contact elevation grid calculated by the PlaneTrendSuf program.Figure SEQ Figure \* ARABIC 2: Appearance of the top contact elevation grid calculated by the PlaneTrendSuf program.centercenter00-12706791325Figure 3: Raster calculator setup for generating the "top_resid" residuals grid.Figure SEQ Figure \* ARABIC 3: Raster calculator setup for generating the "top_resid" residuals grid.centercenter050805780405Figure 4: Appearance of the "top_resid" grid.Figure SEQ Figure \* ARABIC 4: Appearance of the "top_resid" grid.centercenter003613157284085Figure 5: Window dialog for the "Reclassify" procedure.Figure SEQ Figure \* ARABIC 5: Window dialog for the "Reclassify" procedure.centercenter0050805780405Figure 6: Appearance of the "Top_Reclass" map.Figure SEQ Figure \* ARABIC 6: Appearance of the "Top_Reclass" map.centercenter0050805780405Figure 7: Results of adding the "top_resid" and "bot_resid" to produce ""Litho_raster".Figure SEQ Figure \* ARABIC 7: Results of adding the "top_resid" and "bot_resid" to produce ""Litho_raster".centercenter00-25405715635Figure 8: SAGA raster contouring procedure window dialog.Figure 8: SAGA raster contouring procedure window dialog.centercenter03003555408295Figure 9: SAGA line simplification window dialog.Figure 9: SAGA line simplification window dialog.centercenter0005890895Figure 10: Creating the starting lithology polygon.Figure 10: Creating the starting lithology polygon.centercenter0005890895Figure 11: Sliced “Lithology" polygons.Figure 11: Sliced “Lithology" polygons.centercenter00Figure 12: Style settings for the "Lithology" polygon layer.-6355735320Figure 13: Calculation query for the "area" field of the "Lithology" polygons.Figure SEQ Figure \* ARABIC 13: Calculation query for the "area" field of the "Lithology" polygons.centercenter0centercenter0-31756297930Figure SEQ Figure \* ARABIC 14: Example problem map composition layout.Figure SEQ Figure \* ARABIC 14: Example problem map composition layout. ................
................

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

Google Online Preview   Download