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
Locked
