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

[gene]

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)

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)

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

6) result of interpolation of the surface (limit formations A-B) with r.surf.nnbathy (l algorithm)

Hoping that it may be useful to someone. The only (small) problem is you get lots of layers.

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