GRASS GSoC 2013 An Implementation of Horizon Based Stratigraphy for GRASS

(See also other GRASS GSoC 2013 projects)

Student Name:	Tim Bailey, Humboldt State University
Organization:	OSGeo - Open Source Geospatial Foundation
Mentor Name:	Mentor: Benjamin Ducke Backup Mentor: Soeren Gebbert
Title:	An Implementation of Horizon Based Stratigraphy

Abstract

Horizon descriptions are a ubiquitous geographic abstraction for describing three dimensional phenomena in the earth sciences. One of the major trends in modern computational geography is to render continuous fields as discrete finite elements. Substantial progress has been made, applying two dimensional matrices or rasters, to continuous fields. Essentially horizon descriptions are one dimensional sections of three dimensional bodies, with attribute data assigned to intervals within the section. Applying horizon descriptions to populate discrete three dimensional volumes is an evolutionary extension of these current trends. Different scientific disciplines are drawn to different sets of assumptions about the geographic abstraction, based on the particulars of the phenomenon that they study. Rendering three dimensional objects introduces novel issues. Thus, in building interdisciplinary toolsets like GRASS, it is important to make the assumptions used as transparent as possible. In addition it is important to evaluate the applicability of the method for the particular end use.

Weekly Progress reports

Week 1

During this first week of coding I have been working on input specifications for the horizon descriptions interpolation module. Several early technical challenges emerged. The process for importing horzons turns out to be a bit confounding. Generating input from horizon databases turns out to be non-trivial task. I have done an informal survey of every geologist and soil scientist that I know concerning how they store horizon datasets. In addition I have inspected the storage specifications for several public soils database formats. In parallel to the programming project I have been building three sandboxes to test the methods developed in this project. One is a site level wetland delineation. Another is a database of geotechnical borehole logs that Bob Moskovitz of the California Geologic Survey provided. The third is a random Natural Resources Conservation Service survey subset. My work product next week is going to focus on the interpolation module. The abstract workflow for the module as whole is as follows

Process for generation of r3 voxel grid from a population of xy located 1 dimensional horizon descriptions

1. Import of database horizon descriptions
2. Generation of line vector representing the path of the horizon description in xyz
3. Segmentation of the line vector into n number xyz points
4. query of attribute values to points from attribute database
5. Generation of r3 region
6. Interpolation of point attribute values through geostatistical and/or logical operators onto voxel grid locations

After extensive consultation with my programming mentor Ben Ducke, I am now moving on to the interpolation of vector attributes to voxel cells. One point that Ben articulated was that the interpolation assignment needs to be capable of operating on integer data so that it can work on

Week 2

Hello everyone, this week I am checking in a bit early because I am heading to the California Forest Soils Council meeting in a matter of hours. This week I developed a workspace based on prior wetland mapping of a tidally influenced salt marsh in the Northeast corner of Humboldt Bay, in Northwest California. I went through each of the r3 modules to examine the current capabilities. In addition I have been continuing to develop workspaces, focused on subsurface geologic structures in Simi Valley, and Mountain View California. One of the issues that is emerging are the irregularities in how the borehole logs were recorded. I looked into ordinal logistic regression, as well as Tprogs for handling attribute data. Several of the published Ore body assays, quantified z anisotropy compared to xy. In addition I have continued to work on a generalized solution for subsurface structures from p wave first arrivals on seismic surveys. Next week my primary task is going to be to settle and code on the voxel assignment operator. Finally I have come up with a bit of branding. Since r.horizon is already being used by a viewshed module I propose that my module to be named r3.strata. Kind of nice ring to it.

Week 3

This week I worked on the voxel assignment operator for my module.  I also explored the possibility of using a priori assigned anisotropy with a kernal operator to populate the voxel space. One concern that I am not sure how to manage yet are horizon instances that are only identified in a single sample. Next week I will continue working on these same areas. 

Week 4

This was an interesting week. After thorough reconsideration of the foundational documents for this project as well as reviewing about a dozen other three d interpolation projects, I have settled on a methodology that I believe resolves many of the outstanding issues. The vertical to horizontal anisotropy problems are a rabbit hole that I do not need to resolve during this process. Rather than get into the weeds of benefits and liabilities of different interpolation operators, I think that the best way to move forward is by deterministic flood population of voxels. Clearly different disciplines have different methods and needs for spatial approximations. In my estimation, the end users decisions about what the appropriate geostatistical or logical operator to apply should take place between generation of conforming horizons, and the operation that generates the r3 map. The most important point of this is the population of the r3 voxels will proceed through a deterministic operator, which can be accomplished with existing tools. The various approximations such as geostatistical operations can be accomplished before the r3 map is initiated. The r3 map can essentially be flood filled from points by selecting all the points in a tier and applying v.voronoi, followed by v.to.rast. A data artifact from this process is an arbitrary step pattern caused by sparse data. If you want to smooth transitions you can introduce intermediate horizon approximations based on the logical model of the operators choosing. I have a proof of concept dataset which I manually shepherded through what is at this point a ten step transformation between a conforming horizon descriptions, and an r3 map. Next week I plan on continuing working on the generic wrapper that manages the transformation of horizon intervals to voxels. Also I will get my project wiki page up on very soon.

Week 5

This week I worked on the voronoi operator. Also I have been working on a few vector pre-processing procedures, that will end up as options for the module. One of them extends horizon interval vectors to the limit of the computational region to limit the influence of the highest and lowest data. Next week I will continue to work on the voronoi operator and put examples and the workflows up on the wiki. One thing that has been blocking me has been the visualization tools for volumes. m.nviz.image seems to work fine but the 3d view in the Map display is not working correctly. I did an svn update to see if that improved things.

Week 6

This week I took a slight detour from last weeks plan and worked on implementing strike and dip notation for vector points. Also I worked on the wiki, as well as the setting up the repository for my project. I will add materials to each of these as I edit them through the weekend. Several things are blocking me. I have successfully output r3 maps with appropriate attribute values, according to r3.info however I still have been unable to render them in the map display. Also I still do not have a functioning voronoi operator to manage zones of intermediate influence. Next week I will be back to working on this problem as well as integrating all of my work thus far. Also I have not figured out how to use a surface model as an upper boundary.

Tim

Manual Notes

NAME

r3.strata – Creates a voxel map of stratigraphic data from a point layer that includes strike, dip, x, y, z coordinates, and an attribute column

KEYWORDS

raster, profile, horizons, raster3

SYNOPSIS

r3.strata r3.strata help r3.strata name=string strike= dip= easting= northing= elevation= type= value= xyrange= surface_mask= [--overwrite] [--verbose] [quiet]

FLAGS:

--overwrite --verbose --quiet

PARAMETERS:

• name=string
  • Name of r3 raster to be created
• strike=float
  • Azimuth of the direction of maximum slope.
  • Default: 0.0
• dip=float
  • Default: 0.0
• easting= point column with x coordinate
• northing= point column with y coordinate
• elevation= point column with z coordinate
• type=string
  • Type of raster to be created.
  • Options: int, float, double
• Value= The point column with attribute
• One of the specific objectives is to provide
• xyrange= parameter limits influence of

surface_mask= A conventional digital elevation model that limits the output of the three dimensional interpolator to subsurface horizons.

EXAMPLE

r3.strata

AUTHOR

Tim Bailey, Humboldt State University based on r.plane and r3.neighbors