For those interested, following my question of "automate v.extrude from a table ?" (<a href="http://osgeo-org.1560.n6.nabble.com/automate-v-extrude-from-a-table-td4978722.html">http://osgeo-org.1560.n6.nabble.com/automate-v-extrude-from-a-table-td4978722.html</a>), I found a solution in pure Python (consulting <a href="href=%22http://code.google.com/p/postgis-grass-r-py/wiki/0003_01_PythonForGrassGis%22%3Ehttp://code.google.com/p/postgis-grass-r-py/wiki/0003_01_PythonForGrassGis">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</a>). <br>
<br><b><span id="result_box" class="" lang="en"><span class="hps">My problem is</span><span class=""> to represent the</span> <span class="hps">3D</span> <span class="hps">boreholes</span> <span class="hps">in</span> <span class="hps">nviz</span></span> <span id="result_box" class="" lang="en"><span class="hps">without </span></span><span id="result_box" class="" lang="en"> <span class="hps">manually </span></span><span id="result_box" class="" lang="en"><span class="hps">repeating</span> <span class="hps">the v.extrude command </span><span class="hps">n times</span></span>.</b><br>
<br>My solution in Python (<span id="result_box" class="" lang="en"><span class="hps">that works for me)</span></span>:<br> <br>1) I create a buffer around each point locating a borehole to obtain an area layer (circular)<br>
<br><div style="text-align:center"><img alt="Images intégrées 4" src="cid:ii_137b242987905f46" height="134" width="200"><br></div><br>2) I process this layer in Python (It is also possible to run it in the Python shell of the Layer Manager)<br>
<br><div style="margin-left:80px"> <b># import modules (for the pure Python script</b>)<br> import sys, os, numpy, argparse<br> import grass.script as g<br> import grass.script.setup as gsetup <br> <b># init</b><br>
gisbase = os.environ['GISBASE']<br> gisdb="/Users/me/grassdata"<br> location="geol"<br> mapset="mymapset"<br> gsetup.init(gisbase, gisdb, location, mapset)<br><br> <b># read the depths of the limits of two formations (A and B) in the boreholes</b> <b>table</b><br>
inmap = "boreholes" # the area layer <br> val = g.read_command('v.db.select', map = inmap, columns = "A,B",fs=',', flags = 'c')<br> levels=val.split('\n')<br> levels=[i.split(',') for i in levels]<br>
print levels<br> <i>[['890', '850'], ['1040', '830'],....]</i><br> <b># read the identification code of the boreholes</b> <b>in the table</b><br> test= g.read_command('v.db.select', map = inmap, columns = "IDENT",fs=',', flags = 'c')<br>
ident = test.split('\n')<br> ident=[i.split(',') for i in ident]<br> print ident<br> <i> [['3930148'], ['3930243'],.....]</i> <br> <b> # read info from the vector layer (number of boreholes in the layer)</b><br>
info = g.parse_command('<a href="http://v.info">v.info</a>',flags='t', map=inmap,quiet=True)<br> n = int(info['areas'])<br><br> <b># extrude each individual boreholes (one 3D layer for each borehole)</b><br>
for i in range(n-1):<br> name = "bore"+str(ident[i]) # name of the 3D layer <br> # selection/extraction of one element of the table<br> g.run_command('v.extract',input=inmap, list=i+1, output='tmp',type='area', quiet=True, overwrite=True)<br>
# 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)<br> 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)<br>
g.run_command('g.remove', vect='tmp', quiet=True)<br></div><br><b>3) the resulting 3D layers (cylinders) in nviz (thickness of the interval A-B) </b><br><br><div style="text-align:center"><img alt="Images intégrées 1" src="cid:ii_137b2182d1aef627" height="156" width="420"><br>
</div><br>4) <b>it is not possible to create a single layer because v.patch writes always 2D layers. For example:</b><br> <br><div style="text-align:left;margin-left:80px"> <b># creation of a layer</b><br>
g.run_command('v.edit', tool='create', map=outmap)<br> <br> <b># extrude each borehole and add to the layer (v.patch)</b><br> for i in range(n-1):<br> g.run_command('v.extract',input=inmap, list=i+1, output='tmp',type='area', quiet=True, overwrite=True)<br>
g.run_command('v.extrude', input='tmp', <b>output='tmp_extr'</b>, zshift=float(levels[i][1]),height= float(levelst[i][0])-float(levels[i][1]),overwrite=True)<br> <b>g.run_command('v.patch', flags='a', input='tmp_extr', output=outmap, quiet=True, overwrite=True)</b><br>
g.run_command('g.remove', vect='tmp,tmp_extr', quiet=True)<br></div><br>gives the result "<b>WARNING: The output map is not 3D</b>"<br><br>5) <span id="result_box" class="" lang="en"><span class="hps">I repeat</span> <span class="hps">the process</span> <span class="hps">for the other</span> <span class="hps">intervals</span></span><div style="text-align:center">
<img alt="Images intégrées 2" src="cid:ii_137b2318e8bd5546" height="152" width="420"><br><div style="text-align:left"><span id="result_box" class="" lang="en"><span class="hps">6) result of interpolation of the surface (limit formations A-B) with r.surf.nnbathy (l algorithm)<b><br>
</b></span></span><span id="result_box" class="" lang="en"><span class="hps"><br><br></span></span><div style="text-align:center"><div style="text-align:center"><span id="result_box" class="" lang="en"><span class="hps"><img alt="Images intégrées 5" src="cid:ii_137b262c548e1312" height="132" width="420"></span></span><br>
<br><span id="result_box" class="" lang="en"><span class="hps"><br></span></span><div style="text-align:left"><span id="result_box" class="" lang="en"><span class="hps">
Hoping that</span> <span class="hps">it may</span> <span class="hps">be useful to someone. </span></span><span id="result_box" class="" lang="en"><span class="hps">The only</span> <span class="hps atn">(</span><span class="">small)</span> <span class="hps">problem</span> <span class="hps">is you get</span> <span class="hps">lots of</span> <span class="hps">layers.</span></span><br>
</div></div><span id="result_box" class="" lang="en"><span class="hps"></span></span></div></div></div>