3D geological volume modeling (raster 3D): is it really possible ?

[gene]

Hello

Since I do not have Gocad (or others 3D geological modeling softwares), I try to do 3D modeling in geology with GRASS GIS. It works very well for surfaces

• :faults surfaces from contours ( from boreholes and coal mining data):

• Thrust faults from boreholes data and traces of the faults on geological map:

• geological layers and fault (in red) from  boundaries of formations on a geological map (using Eric de Kemp methods in: "Three-dimensional projection of curvilinear geological features through direction cosine interpolation of structural field observations", Computers & Geosciences, Volume 24, n° 3 (1998), pages 269-284 (http://www.sciencedirect.com/science/article/pii/S0098300497000666) or his thesis "THREE-DIMENSIONAL INTEGRATION AND VISUALIZATION OF STRUCTURAL FIELD DATA: TOOLS FOR REGIONAL SUBSURFACE MAPPING" (2000) available at http://bibvir.uqac.ca/theses/12127333/12127333.pdf  and Maciej Tomaszczyk work at http://opengeology.pl/index.php?act=category&category_id=5), see http://www.portailsig.org/content/grass-gis-et-paraview-la-modelisation-geologique-3d-pour-les-demunis-ou-les-enthousiastes (in french)

• interpolation of planar geological layers (v.surf.rst) from boreholes data 

• and cross sections from the model with nviz "cutting planes"

But my problem is to build volumes (raster 3D) from these surfaces.  I read everything that was available on the series of r3.x modules and VTK with Paraview and the work of Sören Gebbert but I do not get any results really satisfactory.

For example:

from to

Is this really possible ?  (it needs lots and lots of memory)

Many thanks

_______________________________________________
grass-user mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-user

[Michael Barton]

It is indeed possible. The 2 important issues are the way you construct the voxels and the way you visualize them.

If you want to have solid layers with sharp upper and lower boundaries (e.g., geologic strata), use r.to.rast3 or r.to.rast3elev. If you want to represent a property that varies continuously vertically (e.g., soil OM or groundwater), you need to create 3D points and do a 3D interpolation. 

Once you have a volume, it just looks like a block unless you specify colors and visualize it in NVIZ with isosurfaces and/or slices. However, even if it "looks" like a block, you can still query it as a 3D volume to extract a geologic layer, for example based on its z coordinate or voxel value.

Michael

On Jun 9, 2012, at 6:04 AM, <[hidden email]>

 wrote:

For example:

from <plangeol.jpg>to <volumo.jpg>

Is this really possible ?  (it needs lots and lots of memory)

Many thanks 

_____________________

C. Michael Barton

Visiting Scientist, Integrated Science Program

National Center for Atmospheric Research &

University Consortium for Atmospheric Research

303-497-2889 (voice)

Director, Center for Social Dynamics & Complexity 

Professor of Anthropology, School of Human Evolution & Social Change

Arizona State University

www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

_______________________________________________
grass-user mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-user

[gene]

I think that the problem is my understanding of  r.to.rast3elev

The problem


I define the 3D resolution (the lower limit of the formation is 540m and the upper 1330m)
   g.region res=200 res3=200  t=1400  b=500 tbres=20
(2d/3d resolution of 200m, top at 1400m, bottom at 500m and a top-bottom-resolution of 20m.

I create two support maps
   r.mapcalc 'one= 1'
   r.mapcalc 'two= 2'

and

r.to.rast3elev --o -l input=one,two elevation=bxltop,bxlbase output=bxl
If I export to vtk the result (with the boreholes) is (with the threshold filter = 1)


The upper surface is correct, the lower not


and and I do not understand why

[gene]

Ok, it works now

with the two surfaces


the volume between the surfaces:


with the base surface:

and the top surface:


But I can not visualize it in nviz:
- I tried with isosurfaces attributes 1 and 2 or 2 only, like the threshold values in Paraview (from  r.mapcalc 'one= 1' and r.mapcalc 'two= 2' )
GRASS 6.4.2 (geol):~ > r3.univar volume20
 100%
total null and non-null cells: 2410200
total null cells: 1703698
Of the non-null cells:
----------------------
n: 706502
minimum: 2
maximum: 2
range: 0
mean: 2
mean of absolute values: 2
standard deviation: 0
variance: 0
variation coefficient: 0 %
sum: 1413004

- I tried with different values ​​associated with surfaces :

GRASS 6.4.2 (geol):~ > r.univar bxlbase
 100%
total null and non-null cells: 48204
total null cells: 8133
Of the non-null cells:
----------------------
n: 40071
minimum: 463.652
maximum: 749.639
range: 285.987
mean: 577.433
mean of absolute values: 577.433
standard deviation: 87.919
variance: 7729.75
variation coefficient: 15.2258 %
sum: 23138333.880065918
GRASS 6.4.2 (geol):~ > r.univar bxlbase
 100%
total null and non-null cells: 48204
total null cells: 8133
Of the non-null cells:
----------------------
n: 40071
minimum: 544.634
maximum: 1081.48
range: 536.846
mean: 842.635
mean of absolute values: 842.635
standard deviation: 126.068
variance: 15893.2
variation coefficient: 14.9612 %
sum: 33765222.471862793

But my nviz remains hopelessly empty.... (only the box outline)