[GRASS-user] v.vol.rst basic questions

GRASS GRASS user
1

Christian Ferreira wrote:

PS: I'm stuck on this since some days, so I would really appreciate
some enlightning from my GRASS colleagues at least for my first
question.

I have been thinking about it, maybe some ideas below...

I was working with volume interpolation with nearly
no problem following the GRASS documentation and the
grassbook for a UTM location.

But working with a Lat-Long location is creating many
doubts in my mind.

stepping back a bit from the 3D, does any one have a comment about
the appropriateness of using v.surf.rst in lat/lon locations?

First, I don't know v.vol.rst really works with
Lat-Long, since all the examples I saw were in UTM
location?

did you get a result? (be a "try it and find out" experimentalist :slight_smile:

If yes, my location is:

RASS 6.4.0RC4 (Transect11):~ > g.region -p
projection: 3 (Latitude-Longitude)
zone: 0

datum: wgs84
ellipsoid: wgs84
north: 10:57:00.288311S
south: 11:03S
west: 78:49:12W
east: 77:45:13.757688W
nsres: 0:00:00.988219
ewres: 0:00:00.988219

so res ~ 30m?

rows: 364
cols: 3884
cells: 1413776

it is approx 1 degree wide, and even though it is on the edge of UTM 18
you probably can get away with reprojecting half way into the next zone
with acceptable distortion. else just set up a custom transverse mercator
with a lon_0 of 78W instead of 75W to work in.

as you have already guessed, giving x,y,z all using the same units
removes a lot of problems.

And I was setting t=0 b=-2200 tbres=20... to have 100 depths.
But I don't unterstand how should be res3, since I'm not using a
UTM location.
Should I convert, let's say "10 meters" to decimal degrees, and
using the result to set "dmin"?

I don't know, I expect it doesn't know about lat/lon units so it may
consider depths to be in degrees too. but maybe that's ok -- once it
is stretched back to meters the effect will be to induce a barrier to
strong vertical mixing, versus the relative ease of horizontal mixing,
which is in fact somewhat realistic. (at least it's wrong in the right
direction)

NVIZ does know to scale lat/lon maps back to meters for 2.5D raster
maps*, probably it scales 3D volumes too?

[*] that still requires some heavy refinement though

Third, at v.vol.rst, I don't understand how to set "dmin"? Is the
same thing as above?

I don't know, but I'd guess it is "dumb". Also I don't know if dmin
is considered in the voxel cube or just in horizontal slices?

Something you might try is to divide the depths by 1852*60 (approx
meters in a deg lat) to convert the depths into degree units, but
really reprojecting into UTM or so would be my first try. you are
somewhat lucky that you are working near the equator.

Basically I have a raster map with multibeam (bathymetry)
along the 11 degree South parallel, and many CTDs measurements
in the water column.

[for the benefit of the list, a CTD is an instrument that measures
physical properties of the water in a vertical profile lowered
directly below the boat to the sea floor. an atmospheric analogy is
driving cross-country and stopping every hour to deploy a radiosonde
with a retrieval line. it's a problem of high vertical resolution but
very poor horizontal resolution (.5m vs 5000m)]

Please see here to understand what I mean:
http://picasaweb.google.com/chris.for.lists/Map#5333477845413068514

actually the thing that I worry about there is not the x,y with z using
different units (easily fixed by reprojecting to a planimetric
projection), it is that all of your data points are essentially in a
single straight line, and what effect that has on the 3D algorithm?

At the end I would like to create a slice (not exactly a
volume) to show the CTD data over this 3D surface/transect
using nviz.

perhaps remake your axes as distance-along-transect vs. depth? then
treat distance as the x axis and depth as the y axis and make a simple
2D interpolation? (classic matlab contourf(), shading flat or masked
v.surf.rst) then somehow rotate it to be a 1 cell wide 3D raster* to
throw up?

[*] I haven't (yet) written a r3.{in|out}.mat but I suppose it would
not be too hard for you to write a little script to create a single
slice to import with r3.in.ascii (see the help page examples) to rotate
the 2D slice to a vertical one?

check out Helena's Chesapeake Bay interpolations:
http://grass.osgeo.org/screenshots/viz.php (near the end)
http://skagit.meas.ncsu.edu/~helena/gmslab/viz/ches.html

and also maybe some r3.out.vtk stuff with Paraview or Vis5D?

Hamish

2

Hi Hamish,

stepping back a bit from the 3D, does any one have a comment about
the appropriateness of using v.surf.rst in lat/lon locations?

So, I have multibeam data covering the coast of Ecuador and Peru.
And Since this data was in Lat-Long, as well the stations lists, I just tried
to use the same format.

Personally, I try to use Lat-Long because my “working areas” are normally
big.

did you get a result? (be a “try it and find out” experimentalist :slight_smile:

Apparently I have results see bellow, and I see it on nviz.

http://picasaweb.google.com/chris.for.lists/Map#5334573425418213906

±---------------------------------------------------------------------------+
| Layer: cdt_volume2@PERMANENT Date: Mon May 11 16:06:29 2009 |
| Mapset: PERMANENT Login of Creator: christian |
| Location: Transect11 |
| DataBase: /Users/christian/Work/M77/GRASS |
| Title: ( cdt_volume2 ) |

Timestamp: none
Type of Map: 3d cell Number of Categories: 0
Data Type: double
Rows: 359
Columns: 3835
Depths: 10
Total Cells: 13767650
Projection: Latitude-Longitude (zone 0)
N: 10:57:00.288311S S: 11:03S Res: 0:00:01.001982
E: 77:45:13.757688W W: 78:49:12W Res: 0:00:01.000845
T: 0 B: -2200 Res: 220
Range of data: min = 2.56193161 max = 12.80138111
Data Source:
Data Description:
generated by v.vol.rst
Comments:
v.vol.rst input=“ctd_3d” cellinp=“bathy” wcolumn=“temperatur” tensio\
n=40. smooth=0.1 segmax=100 npmin=200 dmin=0.01 wmult=1.0 zmult=0.00\
0278 elev=“cdt_volume2”
±---------------------------------------------------------------------------+

Actually, here I have a very strange thing…

If I run g.list rast3d, I can list my interpoled volume.

But if I try to run “r3.info ctd_volume2” in the command line I get a message
telling “ERROR: Raster map <ctd_volume2> not found”, but if I use the GUI for
r3.info is possible to read the file (like above).

The same happens if I try start nviz with:

“nviz bathy volume=ctd_volume2”

So, here is something wrong, because I’m able to manipulate the file just using the GUI…
I don’t unterstand what? Maybe a bug?

so res ~ 30m?

Yes!

it is approx 1 degree wide, and even though it is on the edge of UTM 18
you probably can get away with reprojecting half way into the next zone
with acceptable distortion. else just set up a custom transverse mercator
with a lon_0 of 78W instead of 75W to work in.

as you have already guessed, giving x,y,z all using the same units
removes a lot of problems.

Aham. So, this means that v.vol.rst is not supposed to be used with Lat-Long?

If yes, I will definitely quit this Lat-Long location and go to UTM.

I don’t know, I expect it doesn’t know about lat/lon units so it may
consider depths to be in degrees too. but maybe that’s ok – once it
is stretched back to meters the effect will be to induce a barrier to
strong vertical mixing, versus the relative ease of horizontal mixing,
which is in fact somewhat realistic. (at least it’s wrong in the right
direction)

NVIZ does know to scale lat/lon maps back to meters for 2.5D raster
maps*, probably it scales 3D volumes too?

[*] that still requires some heavy refinement though

I don’t know, but I’d guess it is “dumb”. Also I don’t know if dmin
is considered in the voxel cube or just in horizontal slices?

Something you might try is to divide the depths by 1852*60 (approx
meters in a deg lat) to convert the depths into degree units, but
really reprojecting into UTM or so would be my first try. you are
somewhat lucky that you are working near the equator.

Well, after trying the whole Friday with nearlly no results, I’m quite convinced
that v.vol.rst is not dumb.

And, I was not able to set very well “dmin”.So, the results were strange.

Here is a try adjusting res3 to 0.000278 (30m), zmult to the same number and dmin to 0.001

http://picasaweb.google.com/chris.for.lists/Map#5334573425418213906

actually the thing that I worry about there is not the x,y with z using
different units (easily fixed by reprojecting to a planimetric
projection), it is that all of your data points are essentially in a
single straight line, and what effect that has on the 3D algorithm?

As long I’m capable to produce a 2D slice along the transect/volume
to show the data, the rest would not matter for the moment.

But I understand that interpolation algorithm is having a hard time when trying to
interpolate the cells away from the transect.

This was a test… the best is probably how you are saying below…

perhaps remake your axes as distance-along-transect vs. depth? then
treat distance as the x axis and depth as the y axis and make a simple
2D interpolation? (classic matlab contourf(), shading flat or masked
v.surf.rst) then somehow rotate it to be a 1 cell wide 3D raster* to
throw up?

[*] I haven’t (yet) written a r3.{in|out}.mat but I suppose it would
not be too hard for you to write a little script to create a single
slice to import with r3.in.ascii (see the help page examples) to rotate
the 2D slice to a vertical one?

Well, is not better (or possible) to set in GRASS a temporary region so thin
(in the N-S axis) to create this 2D slice?

Otherwise, I will simply do it like you said before… creating a 1 cell wide 3D raster.

PS: I using for this tests GRASS in a Macbook/Intel on Mac OS X 10.5, and using
the binaries from http://www.kyngchaos.com/software:grass/

check out Helena’s Chesapeake Bay interpolations:
http://grass.osgeo.org/screenshots/viz.php (near the end)
http://skagit.meas.ncsu.edu/~helena/gmslab/viz/ches.html

and also maybe some r3.out.vtk stuff with Paraview or Vis5D?

I’m trying to simply do all using GRASS… maybe is not the best choice.

Regards,
Chris

–
Christian Ferreira
Oceanographer
IFM-GEOMAR Leibniz-Institute of Marine Sciences
Kiel - Germany

Poseidon Linux team
http://www.poseidonlinux.org