[GRASS-user] Generate 3D raster from two surfaces

GRASS GRASS user
1

I want to generate a 3D raster (and display in nviz) that is the volume
between two 2D elevation rasters. The application is to visualize the
subsuraface area that will be subject to environmental remediation.

The values of the 3D raster will be trivial. They will either be "1"
(or anything) or null.

The 2D rasters have elevation as their cell values.

Can r3.mapcalc be used for this? Can I get grass to recognize the
raster value of the 2D raster as the z coordinate? If I could do that,
it would be as simple as assigning the new 3D raster a value of 1 so
long as x and y were not null, and that z1 > z2.

Or is there some other easier way to do this?

Thanks.

J.S.

2

J.S. wrote:

I want to generate a 3D raster (and display in nviz) that is
the volume between two 2D elevation rasters. The application
is to visualize the subsuraface area that will be subject to
environmental remediation.

maybe this summary helps,
  http://grass.osgeo.org/wiki/Help_with_3D

Hamish

3

Hamish:

I looked at that before and the only thing I could think of was to use
r.to.rast3 by first generating a whole series of horizontal surfaces
with r.mapcalc, one for each foot of elevation, with each cell getting a
value of "1" if the "elevation" for that particular new surface lied
between the corresponding cell values of the top and bottom surfaces.
Then taking those whole bunch of surfaces and feeding them to
r.to.rast3. But that seemed pretty inelegant. But it may be only way
to do it. I would only be about 25 rasters to generate and then merge.

J.S.

On Sun, 20 Feb 2011 17:14:39 -0800 (PST), "Hamish" <hamish_b@yahoo.com>
said:

J.S. wrote:
> I want to generate a 3D raster (and display in nviz) that is
> the volume between two 2D elevation rasters. The application
> is to visualize the subsuraface area that will be subject to
> environmental remediation.

maybe this summary helps,
  http://grass.osgeo.org/wiki/Help_with_3D

Hamish

4

Hi,
You can try r.to.rast3elev which will do exactly what you need.

Best
Soeren

Am 21.02.2011 02:38 schrieb <grass@sundquist.imapmail.org>:

Hamish:

I looked at that before and the only thing I could think of was to use
r.to.rast3 by first generating a whole series of horizontal surfaces
with r.mapcalc, one for each foot of elevation, with each cell getting a
value of “1” if the “elevation” for that particular new surface lied
between the corresponding cell values of the top and bottom surfaces.
Then taking those whole bunch of surfaces and feeding them to
r.to.rast3. But that seemed pretty inelegant. But it may be only way
to do it. I would only be about 25 rasters to generate and then merge.

J.S.

On Sun, 20 Feb 2011 17:14:39 -0800 (PST), “Hamish” <hamish_b@yahoo.com>
said:

J.S. wrote:

I want to generate a 3D raster (and display in nviz) that is
the volume between two 2D elevation rasters. The application
is to visualize the subsuraface area that will be subject to
environmental remediation.

maybe this summary helps,
http://grass.osgeo.org/wiki/Help_with_3D

Hamish

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

5

I admit I am not a pro (but I have devoured both the 2nd and 3rd
editions of "the book", and have been using GRASS on and of for about 5
years now), but I am lost on this command.

Here is part of what I am trying to do: Let's try showing just a 3D
raster of the space between the ground surface and the bottom (I have a
few otehr more intersting applications, too, but this is the simplest).

I put GTiff exports of these files here:
http://public.sundquist.imapmail.org/bottom_2.tif and
http://public.sundquist.imapmail.org/ground_surface.tif

I want to generate a 3D raster showing the space between these two
surfaces.

I tried running this command:

r.to.rast3elev input=ground_surface,AREA_DSM
elevation=ground_surface,AREA_DSM output=3d1

and got no errors but got essentially nothing:

r3.univar input=3d1@PERMANENT
total null and non-null cells: 546
total null cells: 469
Of the non-null cells:
----------------------
n: 77
minimum: 1
maximum: 1
range: 0
mean: 1
mean of absolute values: 1
standard deviation: 0
variance: 0
variation coefficient: 0 %
sum: 77

I added to input files because when I tried just one, it said the number
of input and elevation files must match. Note, I do not care about the
actual values in the 3D raster, just whether they are null or not.

Sorry if this is a basic question, but if I can get this generated, it
would be very helpful for explaining something to a client.

Thanks.

J.S.

On Mon, 21 Feb 2011 06:06:37 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hi,
You can try r.to.rast3elev which will do exactly what you need.

Best
Soeren

Am 21.02.2011 02:38 schrieb <grass@sundquist.imapmail.org>:
> Hamish:
>
> I looked at that before and the only thing I could think of was to use
> r.to.rast3 by first generating a whole series of horizontal surfaces
> with r.mapcalc, one for each foot of elevation, with each cell getting a
> value of "1" if the "elevation" for that particular new surface lied
> between the corresponding cell values of the top and bottom surfaces.
> Then taking those whole bunch of surfaces and feeding them to
> r.to.rast3. But that seemed pretty inelegant. But it may be only way
> to do it. I would only be about 25 rasters to generate and then merge.
>
> J.S.
>
> On Sun, 20 Feb 2011 17:14:39 -0800 (PST), "Hamish" <hamish_b@yahoo.com>
> said:
>> J.S. wrote:
>> > I want to generate a 3D raster (and display in nviz) that is
>> > the volume between two 2D elevation rasters. The application
>> > is to visualize the subsuraface area that will be subject to
>> > environmental remediation.
>>
>> maybe this summary helps,
>> http://grass.osgeo.org/wiki/Help_with_3D
>>
>>
>> Hamish
>>
>>
>>
>>
> _______________________________________________
> grass-user mailing list
> grass-user@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user

6

Hi,
make sure your 3d region is set correctly, like:

g.region res=200 res3=200 t=2000 b=0 tbres=20

This will set a region with a 2d/3d resolution of 200m, the top at
2000m, the bottom at 0m and a top-bottom-resolution of 20m.

In this case you will have 100 layer.

Now fill the voxel cells below your surface and bottom with different
values, first we
create two helper maps which represent the values below the elevation maps.

r.mapcalc "one = 1"
r.mapcalc "two = 2"

Now start the module and set the elevation maps from top to bottom,
the flag -l will fill the cells below the elevation maps with the
assigned values.

r.to.rast3elev -l input=one,two elevation=surface,bottom output=volumemap

Your region of interest in the output volume map should have value 1.

You can use paraview to visualize the voxel map. Use r3.out.vtk to
export the voxel map.
Use the threshold filter in paraview to extract all cells with value 1.

Example:
r3.out.vtk input=volumemap output=/tmp/volumemap.vtk
paraview --data=/tmp/volumemap.vtk

Here are two (a bit old) presentations about 3d data handling with
grass and paraview:
http://www-pool.math.tu-berlin.de/~soeren/grass/files/LausanneConferencePresent_soeren_handouts.pdf
http://www-pool.math.tu-berlin.de/~soeren/grass/files/3DWorkshop_soeren_handouts.pdf

Best regards
Soeren

2011/2/21 <grass@sundquist.imapmail.org>:

I admit I am not a pro (but I have devoured both the 2nd and 3rd
editions of "the book", and have been using GRASS on and of for about 5
years now), but I am lost on this command.

Here is part of what I am trying to do: Let's try showing just a 3D
raster of the space between the ground surface and the bottom (I have a
few otehr more intersting applications, too, but this is the simplest).

I put GTiff exports of these files here:
http://public.sundquist.imapmail.org/bottom_2.tif and
http://public.sundquist.imapmail.org/ground_surface.tif

I want to generate a 3D raster showing the space between these two
surfaces.

I tried running this command:

r.to.rast3elev input=ground_surface,AREA_DSM
elevation=ground_surface,AREA_DSM output=3d1

and got no errors but got essentially nothing:

r3.univar input=3d1@PERMANENT
total null and non-null cells: 546
total null cells: 469
Of the non-null cells:
----------------------
n: 77
minimum: 1
maximum: 1
range: 0
mean: 1
mean of absolute values: 1
standard deviation: 0
variance: 0
variation coefficient: 0 %
sum: 77

I added to input files because when I tried just one, it said the number
of input and elevation files must match. Note, I do not care about the
actual values in the 3D raster, just whether they are null or not.

Sorry if this is a basic question, but if I can get this generated, it
would be very helpful for explaining something to a client.

Thanks.

J.S.

On Mon, 21 Feb 2011 06:06:37 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hi,
You can try r.to.rast3elev which will do exactly what you need.

Best
Soeren

Am 21.02.2011 02:38 schrieb <grass@sundquist.imapmail.org>:
> Hamish:
>
> I looked at that before and the only thing I could think of was to use
> r.to.rast3 by first generating a whole series of horizontal surfaces
> with r.mapcalc, one for each foot of elevation, with each cell getting a
> value of "1" if the "elevation" for that particular new surface lied
> between the corresponding cell values of the top and bottom surfaces.
> Then taking those whole bunch of surfaces and feeding them to
> r.to.rast3. But that seemed pretty inelegant. But it may be only way
> to do it. I would only be about 25 rasters to generate and then merge.
>
> J.S.
>
> On Sun, 20 Feb 2011 17:14:39 -0800 (PST), "Hamish" <hamish_b@yahoo.com>
> said:
>> J.S. wrote:
>> > I want to generate a 3D raster (and display in nviz) that is
>> > the volume between two 2D elevation rasters. The application
>> > is to visualize the subsuraface area that will be subject to
>> > environmental remediation.
>>
>> maybe this summary helps,
>> http://grass.osgeo.org/wiki/Help_with_3D
>>
>>
>> Hamish
>>
>>
>>
>>
> _______________________________________________
> grass-user mailing list
> grass-user@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

7

Thanks for the help, but still no luck.

I ran as you suggested:

GRASS 6.4> g.region res=1 res3=1 t=75 b=0 tbres=1
GRASS 6.4> r.mapcalc "one = 1"
GRASS 6.4> r.mapcalc "two = 2"
GRASS 6.4> r.to.rast3elev -l input=one,two elevation=contours,bottom_2
output=volumemap
Creating 3D raster map
100%
100%
GRASS 6.4> r3.out.vtk input=volumemap output=/tmp/volumemap.vtk
100%
GRASS 6.4>

Note, I used just the 1-foot grid (it is a US State Plane coordinate
system) which perhaps is too fine, but I have a relatively small area
(800 x 900 ft) and only need to evaluate 75 feet in thickness and it
still works (just slow). It gave me a 3D raster:

GRASS 6.4> r3.univar input=volumemap
total null and non-null cells: 54000000
total null cells: 40596983

Of the non-null cells:
----------------------
n: 13403017
minimum: 1
maximum: 2
range: 1
mean: 1.51073
mean of absolute values: 1.51073
standard deviation: 0.499885
variance: 0.249885
variation coefficient: 33.089 %
sum: 20248304
100%
----------------------

The overall volume (which would be the number of non-null cells, in
cubic feet) looks a bit high, but may be correct.

However, nothing shows up in paraview or nviz.

Screenshot: http://public.sundquist.imapmail.org/paraview-1.png

It's a pretty by vtk file (~1 GB), but that shouldn't matter.

In paraview, the "Filter" menu is greyed out. I went to tools>create
custom filter but there were no "objects" or "properties" to select. I
am not too familiar with paraview.

Thanks again for the help.

J.S.

On Mon, 21 Feb 2011 15:57:59 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hi,
make sure your 3d region is set correctly, like:

g.region res=200 res3=200 t=2000 b=0 tbres=20

This will set a region with a 2d/3d resolution of 200m, the top at
2000m, the bottom at 0m and a top-bottom-resolution of 20m.

In this case you will have 100 layer.

Now fill the voxel cells below your surface and bottom with different
values, first we
create two helper maps which represent the values below the elevation
maps.

r.mapcalc "one = 1"
r.mapcalc "two = 2"

Now start the module and set the elevation maps from top to bottom,
the flag -l will fill the cells below the elevation maps with the
assigned values.

r.to.rast3elev -l input=one,two elevation=surface,bottom output=volumemap

Your region of interest in the output volume map should have value 1.

You can use paraview to visualize the voxel map. Use r3.out.vtk to
export the voxel map.
Use the threshold filter in paraview to extract all cells with value 1.

Example:
r3.out.vtk input=volumemap output=/tmp/volumemap.vtk
paraview --data=/tmp/volumemap.vtk

Here are two (a bit old) presentations about 3d data handling with
grass and paraview:
http://www-pool.math.tu-berlin.de/~soeren/grass/files/LausanneConferencePresent_soeren_handouts.pdf
http://www-pool.math.tu-berlin.de/~soeren/grass/files/3DWorkshop_soeren_handouts.pdf

Best regards
Soeren

2011/2/21 <grass@sundquist.imapmail.org>:
> I admit I am not a pro (but I have devoured both the 2nd and 3rd
> editions of "the book", and have been using GRASS on and of for about 5
> years now), but I am lost on this command.
>
> Here is part of what I am trying to do: Let's try showing just a 3D
> raster of the space between the ground surface and the bottom (I have a
> few otehr more intersting applications, too, but this is the simplest).
>
> I put GTiff exports of these files here:
> http://public.sundquist.imapmail.org/bottom_2.tif and
> http://public.sundquist.imapmail.org/ground_surface.tif
>
> I want to generate a 3D raster showing the space between these two
> surfaces.
>
> I tried running this command:
>
> r.to.rast3elev input=ground_surface,AREA_DSM
> elevation=ground_surface,AREA_DSM output=3d1
>
> and got no errors but got essentially nothing:
>
> r3.univar input=3d1@PERMANENT
> total null and non-null cells: 546
> total null cells: 469
> Of the non-null cells:
> ----------------------
> n: 77
> minimum: 1
> maximum: 1
> range: 0
> mean: 1
> mean of absolute values: 1
> standard deviation: 0
> variance: 0
> variation coefficient: 0 %
> sum: 77
>
> I added to input files because when I tried just one, it said the number
> of input and elevation files must match. Note, I do not care about the
> actual values in the 3D raster, just whether they are null or not.
>
> Sorry if this is a basic question, but if I can get this generated, it
> would be very helpful for explaining something to a client.
>
> Thanks.
>
> J.S.
>
> On Mon, 21 Feb 2011 06:06:37 +0100, "Soeren Gebbert"
> <soerengebbert@googlemail.com> said:
>> Hi,
>> You can try r.to.rast3elev which will do exactly what you need.
>>
>> Best
>> Soeren
>>
>> Am 21.02.2011 02:38 schrieb <grass@sundquist.imapmail.org>:
>> > Hamish:
>> >
>> > I looked at that before and the only thing I could think of was to use
>> > r.to.rast3 by first generating a whole series of horizontal surfaces
>> > with r.mapcalc, one for each foot of elevation, with each cell getting a
>> > value of "1" if the "elevation" for that particular new surface lied
>> > between the corresponding cell values of the top and bottom surfaces.
>> > Then taking those whole bunch of surfaces and feeding them to
>> > r.to.rast3. But that seemed pretty inelegant. But it may be only way
>> > to do it. I would only be about 25 rasters to generate and then merge.
>> >
>> > J.S.
>> >
>> > On Sun, 20 Feb 2011 17:14:39 -0800 (PST), "Hamish" <hamish_b@yahoo.com>
>> > said:
>> >> J.S. wrote:
>> >> > I want to generate a 3D raster (and display in nviz) that is
>> >> > the volume between two 2D elevation rasters. The application
>> >> > is to visualize the subsuraface area that will be subject to
>> >> > environmental remediation.
>> >>
>> >> maybe this summary helps,
>> >> http://grass.osgeo.org/wiki/Help_with_3D
>> >>
>> >>
>> >> Hamish
>> >>
>> >>
>> >>
>> >>
>> > _______________________________________________
>> > grass-user mailing list
>> > grass-user@lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/grass-user
> _______________________________________________
> grass-user mailing list
> grass-user@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
>
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

8

Hi,
try:

r3.out.vtk input=volumemap output=/tmp/volumemap.vtk null=0.0

Open it in paraview and press the Apply button (thats important).

Then choose style -> representation -> surface in the Display tab and
then color by -> volumemap in the color tab. I suggest you use a
recent paraview version?

To compute the number of cells between the surface and the bottom use
r3.mapcalc:

r3.mapcalc "valid_cells = if(volumemap == 1, 1, null())"

r3.univar valid_cells

To avoid large vtk files while testing, you can reduce the region
resolution with g.region.

Best regards
Soeren

2011/2/21 <grass@sundquist.imapmail.org>:

Thanks for the help, but still no luck.

I ran as you suggested:

GRASS 6.4> g.region res=1 res3=1 t=75 b=0 tbres=1
GRASS 6.4> r.mapcalc "one = 1"
GRASS 6.4> r.mapcalc "two = 2"
GRASS 6.4> r.to.rast3elev -l input=one,two elevation=contours,bottom_2
output=volumemap
Creating 3D raster map
100%
100%
GRASS 6.4> r3.out.vtk input=volumemap output=/tmp/volumemap.vtk
100%
GRASS 6.4>

Note, I used just the 1-foot grid (it is a US State Plane coordinate
system) which perhaps is too fine, but I have a relatively small area
(800 x 900 ft) and only need to evaluate 75 feet in thickness and it
still works (just slow). It gave me a 3D raster:

GRASS 6.4> r3.univar input=volumemap
total null and non-null cells: 54000000
total null cells: 40596983

Of the non-null cells:
----------------------
n: 13403017
minimum: 1
maximum: 2
range: 1
mean: 1.51073
mean of absolute values: 1.51073
standard deviation: 0.499885
variance: 0.249885
variation coefficient: 33.089 %
sum: 20248304
100%
----------------------

The overall volume (which would be the number of non-null cells, in
cubic feet) looks a bit high, but may be correct.

However, nothing shows up in paraview or nviz.

Screenshot: http://public.sundquist.imapmail.org/paraview-1.png

It's a pretty by vtk file (~1 GB), but that shouldn't matter.

In paraview, the "Filter" menu is greyed out. I went to tools>create
custom filter but there were no "objects" or "properties" to select. I
am not too familiar with paraview.

Thanks again for the help.

J.S.

On Mon, 21 Feb 2011 15:57:59 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hi,
make sure your 3d region is set correctly, like:

g.region res=200 res3=200 t=2000 b=0 tbres=20

This will set a region with a 2d/3d resolution of 200m, the top at
2000m, the bottom at 0m and a top-bottom-resolution of 20m.

In this case you will have 100 layer.

Now fill the voxel cells below your surface and bottom with different
values, first we
create two helper maps which represent the values below the elevation
maps.

r.mapcalc "one = 1"
r.mapcalc "two = 2"

Now start the module and set the elevation maps from top to bottom,
the flag -l will fill the cells below the elevation maps with the
assigned values.

r.to.rast3elev -l input=one,two elevation=surface,bottom output=volumemap

Your region of interest in the output volume map should have value 1.

You can use paraview to visualize the voxel map. Use r3.out.vtk to
export the voxel map.
Use the threshold filter in paraview to extract all cells with value 1.

Example:
r3.out.vtk input=volumemap output=/tmp/volumemap.vtk
paraview --data=/tmp/volumemap.vtk

Here are two (a bit old) presentations about 3d data handling with
grass and paraview:
http://www-pool.math.tu-berlin.de/~soeren/grass/files/LausanneConferencePresent_soeren_handouts.pdf
http://www-pool.math.tu-berlin.de/~soeren/grass/files/3DWorkshop_soeren_handouts.pdf

Best regards
Soeren

2011/2/21 <grass@sundquist.imapmail.org>:
> I admit I am not a pro (but I have devoured both the 2nd and 3rd
> editions of "the book", and have been using GRASS on and of for about 5
> years now), but I am lost on this command.
>
> Here is part of what I am trying to do: Let's try showing just a 3D
> raster of the space between the ground surface and the bottom (I have a
> few otehr more intersting applications, too, but this is the simplest).
>
> I put GTiff exports of these files here:
> http://public.sundquist.imapmail.org/bottom_2.tif and
> http://public.sundquist.imapmail.org/ground_surface.tif
>
> I want to generate a 3D raster showing the space between these two
> surfaces.
>
> I tried running this command:
>
> r.to.rast3elev input=ground_surface,AREA_DSM
> elevation=ground_surface,AREA_DSM output=3d1
>
> and got no errors but got essentially nothing:
>
> r3.univar input=3d1@PERMANENT
> total null and non-null cells: 546
> total null cells: 469
> Of the non-null cells:
> ----------------------
> n: 77
> minimum: 1
> maximum: 1
> range: 0
> mean: 1
> mean of absolute values: 1
> standard deviation: 0
> variance: 0
> variation coefficient: 0 %
> sum: 77
>
> I added to input files because when I tried just one, it said the number
> of input and elevation files must match. Note, I do not care about the
> actual values in the 3D raster, just whether they are null or not.
>
> Sorry if this is a basic question, but if I can get this generated, it
> would be very helpful for explaining something to a client.
>
> Thanks.
>
> J.S.
>
> On Mon, 21 Feb 2011 06:06:37 +0100, "Soeren Gebbert"
> <soerengebbert@googlemail.com> said:
>> Hi,
>> You can try r.to.rast3elev which will do exactly what you need.
>>
>> Best
>> Soeren
>>
>> Am 21.02.2011 02:38 schrieb <grass@sundquist.imapmail.org>:
>> > Hamish:
>> >
>> > I looked at that before and the only thing I could think of was to use
>> > r.to.rast3 by first generating a whole series of horizontal surfaces
>> > with r.mapcalc, one for each foot of elevation, with each cell getting a
>> > value of "1" if the "elevation" for that particular new surface lied
>> > between the corresponding cell values of the top and bottom surfaces.
>> > Then taking those whole bunch of surfaces and feeding them to
>> > r.to.rast3. But that seemed pretty inelegant. But it may be only way
>> > to do it. I would only be about 25 rasters to generate and then merge.
>> >
>> > J.S.
>> >
>> > On Sun, 20 Feb 2011 17:14:39 -0800 (PST), "Hamish" <hamish_b@yahoo.com>
>> > said:
>> >> J.S. wrote:
>> >> > I want to generate a 3D raster (and display in nviz) that is
>> >> > the volume between two 2D elevation rasters. The application
>> >> > is to visualize the subsuraface area that will be subject to
>> >> > environmental remediation.
>> >>
>> >> maybe this summary helps,
>> >> http://grass.osgeo.org/wiki/Help_with_3D
>> >>
>> >>
>> >> Hamish
>> >>
>> >>
>> >>
>> >>
>> > _______________________________________________
>> > grass-user mailing list
>> > grass-user@lists.osgeo.org
>> > http://lists.osgeo.org/mailman/listinfo/grass-user
> _______________________________________________
> grass-user mailing list
> grass-user@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
>
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

9

Soeren:

Well, that sort-of worked (at least after I increased the cell size so
paraview didn't crash), but it didn't produce exactly what I was
looking for:

http://public.sundquist.imapmail.org/paraview-2.png

It is basically the rectangular block with the 800 x 900 x 75 foot
dimensions, with the plan-view of the treatment area painted on the top
and bottom (gray one side, red the other side).

And actually, the shape shown on each side is not exactly the treatment
area. For this test example, the top and bottom surfaces have the same
footprint (i.e. plan view).

Here is the top: http://public.sundquist.imapmail.org/top.png

and here is the bottom: http://public.sundquist.imapmail.org/bottom.png

The bottom side of the "cube" (red outline, no screen shot) is similar
in shape to the input rasters, allowing for increased cell size.

But the top side (grey, in screenshot above) has some sections missing
(for example, the rest of the NW corner). It shouldn't have had
sectiosn truncated, since my region was set from 0 to 75 feet elevation:

GRASS 6.4> g.region res=1 res3=10 t=75 b=0 tbres=10

and the top and bottom are within these bounds:

GRASS 6.4> r.univar contours
total null and non-null cells: 720000
total null cells: 525599

Of the non-null cells:
----------------------
n: 194401
minimum: 55.2507
maximum: 71.2332
range: 15.9825
mean: 68.4245
mean of absolute values: 68.4245
standard deviation: 2.02751
variance: 4.11078
variation coefficient: 2.96313 %
sum: 13301789.790035248
100%

GRASS 6.4> r.univar bottom_2
total null and non-null cells: 720000
total null cells: 525599

Of the non-null cells:
----------------------
n: 194401
minimum: 15.5567
maximum: 45.1165
range: 29.5598
mean: 34.7136
mean of absolute values: 34.7136
standard deviation: 3.25303
variance: 10.5822
variation coefficient: 9.37105 %
sum: 6748356.0299701691
100%

Does it matter that I had different values for res and res3?

J.S.

On Mon, 21 Feb 2011 18:30:36 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hi, try:

r3.out.vtk input=volumemap output=/tmp/volumemap.vtk null=0.0

Open it in paraview and press the Apply button (thats important).

Then choose style -> representation -> surface in the Display tab and
then color by -> volumemap in the color tab. I suggest you use a
recent paraview version?

To compute the number of cells between the surface and the bottom use
r3.mapcalc:

r3.mapcalc "valid_cells = if(volumemap == 1, 1, null())"

r3.univar valid_cells

To avoid large vtk files while testing, you can reduce the region
resolution with g.region.

Best regards Soeren

10

Hello,
IMHO the result on the top of the volume data in paraview looks
correct, because the top elevation map as lower values as 75ft at the
NW side.
I strongly recommend that you use the threshold filter and the cutting
tools of paraview to have a deeper insight of the result.

In case the 2d resolution is different from the 3d resolution
internally the 3d resolution in x,y plane is used to read the 2d
raster maps.

Best regards
Soeren

2011/2/21 <grass@sundquist.imapmail.org>:

Soeren:

Well, that sort-of worked (at least after I increased the cell size so
paraview didn't crash), but it didn't produce exactly what I was
looking for:

http://public.sundquist.imapmail.org/paraview-2.png

It is basically the rectangular block with the 800 x 900 x 75 foot
dimensions, with the plan-view of the treatment area painted on the top
and bottom (gray one side, red the other side).

And actually, the shape shown on each side is not exactly the treatment
area. For this test example, the top and bottom surfaces have the same
footprint (i.e. plan view).

Here is the top: http://public.sundquist.imapmail.org/top.png

and here is the bottom: http://public.sundquist.imapmail.org/bottom.png

The bottom side of the "cube" (red outline, no screen shot) is similar
in shape to the input rasters, allowing for increased cell size.

But the top side (grey, in screenshot above) has some sections missing
(for example, the rest of the NW corner). It shouldn't have had
sectiosn truncated, since my region was set from 0 to 75 feet elevation:

GRASS 6.4> g.region res=1 res3=10 t=75 b=0 tbres=10

and the top and bottom are within these bounds:

GRASS 6.4> r.univar contours
total null and non-null cells: 720000
total null cells: 525599

Of the non-null cells:
----------------------
n: 194401
minimum: 55.2507
maximum: 71.2332
range: 15.9825
mean: 68.4245
mean of absolute values: 68.4245
standard deviation: 2.02751
variance: 4.11078
variation coefficient: 2.96313 %
sum: 13301789.790035248
100%

GRASS 6.4> r.univar bottom_2
total null and non-null cells: 720000
total null cells: 525599

Of the non-null cells:
----------------------
n: 194401
minimum: 15.5567
maximum: 45.1165
range: 29.5598
mean: 34.7136
mean of absolute values: 34.7136
standard deviation: 3.25303
variance: 10.5822
variation coefficient: 9.37105 %
sum: 6748356.0299701691
100%

Does it matter that I had different values for res and res3?

J.S.

On Mon, 21 Feb 2011 18:30:36 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hi, try:

r3.out.vtk input=volumemap output=/tmp/volumemap.vtk null=0.0

Open it in paraview and press the Apply button (thats important).

Then choose style -> representation -> surface in the Display tab and
then color by -> volumemap in the color tab. I suggest you use a
recent paraview version?

To compute the number of cells between the surface and the bottom use
r3.mapcalc:

r3.mapcalc "valid_cells = if(volumemap == 1, 1, null())"

r3.univar valid_cells

To avoid large vtk files while testing, you can reduce the region
resolution with g.region.

Best regards Soeren

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

11

On Mon, 21 Feb 2011 21:40:44 +0100, "Soeren Gebbert"
<soerengebbert@googlemail.com> said:

Hello,
IMHO the result on the top of the volume data in paraview looks
correct, because the top elevation map as lower values as 75ft at the
NW side.
I strongly recommend that you use the threshold filter and the cutting
tools of paraview to have a deeper insight of the result.

Thanks. Forgot about the "threshold" suggestion. Works now.

I really appreciate your help on this.

J.S.