[GRASS-user] Automatic 3D geological boreholes representation (automate v.extrude from a table ?): my solution in Python

[Martin Laloux]

For those interested, following my question of "automate v.extrude from a
table ?" (
http://osgeo-org.1560.n6.nabble.com/automate-v-extrude-from-a-table-td4978722.html),
I found a solution in pure Python (consulting
href="http://code.google.com/p/postgis-grass-r-py/wiki/0003_01_PythonForGrassGis">http://code.google.com/p/postgis-grass-r-py/wiki/0003_01_PythonForGrassGis).


*My problem is to represent the 3D boreholes in nviz without  manually
repeating the v.extrude command n times.*

My solution in Python (that works for me):

1) I create a buffer around each point locating a borehole to obtain an
area layer (circular)

[image: Images intégrées 4]

2) I process this layer in Python (It is also possible to run it in the
Python shell of the Layer Manager)

  *# import modules (for the pure Python script*)
   import sys, os, numpy, argparse
   import grass.script as g
   import grass.script.setup as gsetup
    *# init*
   gisbase = os.environ['GISBASE']
   gisdb="/Users/me/grassdata"
   location="geol"
   mapset="mymapset"
   gsetup.init(gisbase, gisdb, location, mapset)

   *# read the depths of the limits of two formations (A and B) in the
boreholes* *table*
   inmap = "boreholes" # the area layer
   val = g.read_command('v.db.select', map = inmap, columns = "A,B",fs=',',
flags = 'c')
   levels=val.split('\n')
   levels=[i.split(',') for i in levels]
   print levels
   *[['890', '850'], ['1040', '830'],....]*
    *# read the identification code of the boreholes* *in the table*
    test= g.read_command('v.db.select', map = inmap, columns =
"IDENT",fs=',', flags = 'c')
    ident =  test.split('\n')
    ident=[i.split(',') for i in ident]
    print ident
  * [['3930148'], ['3930243'],.....]*
  * # read info from the vector layer (number of boreholes in the layer)*
   info = g.parse_command('v.info',flags='t', map=inmap,quiet=True)
   n = int(info['areas'])

   *# extrude each individual boreholes (one 3D layer for each borehole)*
   for i in range(n-1):
         name = "bore"+str(ident[i]) # name of the 3D layer
         # selection/extraction of one element of the table
         g.run_command('v.extract',input=inmap, list=i+1,
output='tmp',type='area', quiet=True, overwrite=True)
         # extrusion with base = limit of B formation (levels[i][1] = 850
for the first), height=thickness A-B (levelst[i][0])-float(levelst[i][1]) =
40 for the first)
         g.run_command('v.extrude', input='tmp', output=name,
zshift=float(levels[i][1]),height=
float(levelst[i][0])-float(levelst[i][1]),overwrite=True)
         g.run_command('g.remove', vect='tmp', quiet=True)

*3) the resulting 3D layers (cylinders) in nviz (thickness of the interval
A-B) *

[image: Images intégrées 1]

4) *it is not possible to create a single layer because v.patch writes
always 2D layers. For example:*

         *# creation of a layer*
         g.run_command('v.edit', tool='create', map=outmap)

         *# extrude each borehole and add to the layer (v.patch)*
         for i in range(n-1):
             g.run_command('v.extract',input=inmap, list=i+1,
output='tmp',type='area', quiet=True, overwrite=True)
             g.run_command('v.extrude', input='tmp', *output='tmp_extr'*,
zshift=float(levels[i][1]),height=
float(levelst[i][0])-float(levels[i][1]),overwrite=True)
             *g.run_command('v.patch', flags='a', input='tmp_extr',
output=outmap, quiet=True, overwrite=True)*
             g.run_command('g.remove', vect='tmp,tmp_extr', quiet=True)

gives the result "*WARNING: The output map is not 3D*"

5) I repeat the process for the other intervals
[image: Images intégrées 2]
6) result of interpolation of the surface (limit formations A-B) with
r.surf.nnbathy (l algorithm)*
*

[image: Images intégrées 5]


Hoping that it may be useful to someone. The only (small) problem is
you get lots
of layers.
-------------- section suivante --------------
Une pi?ce jointe HTML a ?t? nettoy?e...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20120603/74e38a8d/attachment-0001.html>
-------------- section suivante --------------
Une pi?ce jointe autre que texte a ?t? nettoy?e...
Nom: result.png
Type: image/png
Taille: 41915 octets
Desc: non disponible
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20120603/74e38a8d/attachment-0004.png>
-------------- section suivante --------------
Une pi?ce jointe autre que texte a ?t? nettoy?e...
Nom: boreholesa.png
Type: image/png
Taille: 8970 octets
Desc: non disponible
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20120603/74e38a8d/attachment-0005.png>
-------------- section suivante --------------
Une pi?ce jointe autre que texte a ?t? nettoy?e...
Nom: boreholes2.png
Type: image/png
Taille: 9756 octets
Desc: non disponible
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20120603/74e38a8d/attachment-0006.png>
-------------- section suivante --------------
Une pi?ce jointe autre que texte a ?t? nettoy?e...
Nom: boreholes.png
Type: image/png
Taille: 9223 octets
Desc: non disponible
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20120603/74e38a8d/attachment-0007.png>