[GRASS-user] 3D volume interpolation help

Hi all, it's been a while :slight_smile:

I have a student who needs to create volumes from a 3D points file.
The ASCII file looks like this:

741665|7187390|205.67|96|1|1|2|100|0|90
741736.92|7187280.39|193.89|90|2|1|2|92|0|90
741736.92|7187280.39|193.39|50|2|1|2|58|0|90
741665|7187390|191.03|96|1|1|2|100|0|90
741681|7187270|190.85|94|1|1|2|100|0|90

there are 181 entries, which correspond to drill core samples, so
actually we have 9 cores (same XY, different Z)

importing as a 3D vector works fine
v.in.ascii -z input=/home/guano/RQD_2.csv output=rqd_drill z=3

now set up the region:
g.region t=210 b=125 tbres=1 res3=1 -p3
projection: 1 (UTM)
zone: 21
datum: ** unknown (default: WGS84) **
ellipsoid: sam69
north: 7187390
south: 7187270
west: 741665
east: 741761
top: 210.00000000
bottom: 125.00000000
nsres: 1
nsres3: 1
ewres: 1
ewres3: 1
tbres: 1
rows: 120
rows3: 120
cols: 96
cols3: 96
depths: 85
cells: 11520
3dcells: 979200

OK, now on to interpolation (int_1 is the 4th col):

v.vol.rst --overwrite input=RQD_2 wcolumn=int_1 segmax=10 npmin=10
dmin=0.05 elev=teste3d
181 records selected from table
WARNING: strip exists with insufficient data
Processing all selected output files will require 46080 bytes of disk space
for temp files

WARNING: There is less than 700 points for interpolation, no segmentation
         is necessary, to run the program faster, set segmax=700 (see
         manual)
Percent complete: 100%
Finished interpolating

The number of points in vector map is 181
The number of points outside of 2D/3D region 0
The number of points used (after reduction) is 181

OK, let's see what happened:
r3.info teste3d
+----------------------------------------------------------------------------+
| Layer: teste3d Date: Wed Mar 31 15:23:49 2010 |
| Mapset: Vertedouro_Itaipu Login of Creator: guano |
| Location: test_1 |
| DataBase: /home/guano/GIS DataBase |
| Title: ( teste3d ) |
| Timestamp: none |
|----------------------------------------------------------------------------|
| |
| Type of Map: 3d cell Number of Categories: 0 |
| Data Type: double |
| Rows: 120 |
| Columns: 96 |
| Depths: 85 |
| Total Cells: 979200 |
| Projection: UTM (zone 21) |
| N: 7187390 S: 7187270 Res: 1 |
| E: 741761 W: 741665 Res: 1 |
| T: 210 B: 125 Res: 1 |
| Range of data: min = 98 max = 98 |
| |
| Data Source: |
| |
| |
| |
| Data Description: |
| generated by v.vol.rst |
| |
| Comments: |
| v.vol.rst input="RQD_2@Vertedouro_Itaipu" wcolumn="int_1" tension=40\ |
| . smooth=0.1 segmax=10 npmin=10 dmin=0.05 wmult=1.0 zmult=1.0 elev="\ |
| teste3d" |
| |
+----------------------------------------------------------------------------+

So, I got a map with min=max=98 (which is the mean value of that column).

OK, what are we doing wrong? I got to admit I'm new to 3d interpolation.

best

Carlos

v.vol.rst --overwrite input=RQD_2@Vertedouro_Itaipu wcolumn=int_1
segmax=1 npmin=1 dmin=0.005 elev=teste3d

--
Prof. Carlos Henrique Grohmann - Geologist D.Sc.
Institute of Geosciences - Univ. of São Paulo, Brazil
http://www.igc.usp.br/pessoais/guano
Linux User #89721
________________
Can’t stop the signal.

On Wed, Mar 31, 2010 at 8:26 PM, Carlos Grohmann
<carlos.grohmann@gmail.com> wrote:

Hi all, it's been a while :slight_smile:

I have a student who needs to create volumes from a 3D points file.
The ASCII file looks like this:

741665|7187390|205.67|96|1|1|2|100|0|90
741736.92|7187280.39|193.89|90|2|1|2|92|0|90
741736.92|7187280.39|193.39|50|2|1|2|58|0|90
741665|7187390|191.03|96|1|1|2|100|0|90
741681|7187270|190.85|94|1|1|2|100|0|90

there are 181 entries, which correspond to drill core samples, so
actually we have 9 cores (same XY, different Z)

importing as a 3D vector works fine
v.in.ascii -z input=/home/guano/RQD_2.csv output=rqd_drill z=3

now set up the region:
g.region t=210 b=125 tbres=1 res3=1 -p3

Note: I have used (in NC location, since I don't have yours):
g.region t=210 b=125 tbres=1 res3=1 -p3 vect=rqd_drill res=1 -a

projection: 1 (UTM)
zone: 21
datum: ** unknown (default: WGS84) **
ellipsoid: sam69
north: 7187390
south: 7187270
west: 741665
east: 741761
top: 210.00000000
bottom: 125.00000000
nsres: 1
nsres3: 1
ewres: 1
ewres3: 1
tbres: 1
rows: 120
rows3: 120
cols: 96
cols3: 96
depths: 85
cells: 11520
3dcells: 979200

I have different cols and rows but that may be ok for a test.

v.db.select rqd_drill column=int_1
int_1
96
90
50
96
94

OK, now on to interpolation (int_1 is the 4th col):

v.vol.rst --overwrite input=RQD_2 wcolumn=int_1 segmax=10 npmin=10
dmin=0.05 elev=teste3d

warning: above the file name was rqd_drill.

181 records selected from table
WARNING: strip exists with insufficient data
Processing all selected output files will require 46080 bytes of disk space
for temp files

...

OK, let's see what happened:
r3.info teste3d

...

| N: 7187390 S: 7187270 Res: 1 |
| E: 741761 W: 741665 Res: 1 |
| T: 210 B: 125 Res: 1 |
| Range of data: min = 98 max = 98 |

I get:

r3.info teste3d | grep -i Range
| Range of data: min = 67.83296967 max = 94.63486481 |
...

So, I got a map with min=max=98 (which is the mean value of that column).

OK, what are we doing wrong? I got to admit I'm new to 3d interpolation.

... could it be as simple as you used the wrong input file?

I used GRASS 6.4.0RC6 here which appears to work fine.

Markus