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 |
Free forum by Nabble | Edit this page |