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

[Hamish]

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 :)


> 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