[GRASS-user] r3.in.xyz

[Moritz Lennert]

Hi François,

On 10/07/18 04:40, Francois Chartier wrote:
> Hi Moritz,
> 
> The idea is a bit novel.  I am working with over 20000 Boreholes (BH) 
> with soil information (silt, sand, gravel, Tills, and everything in 
> between etc.).  These are not distinct geological contacts with 
> recognizable stratigraphies that can be correlated, ex: this top of 
> gravel corresponds to this top of gravel in this BH etc.  The model 
> domain is approximately 40 km by 40 km and extends from ground surface 
> (DEM) to top of bedrock (surface from geological survey), and thickness 
> varies from approx. 15 m to 100 m.
> As it would take years to correlate all the soil contact elevations, I 
> am approaching the problem from a different angle.  The dataset consists 
> of XYZ and Property (top of soil deposit) at each borehole location 
> (XY).  Along the vertical axis of the borehole, there are multiple top 
> of soil deposits overlying on top of each other depending on the depth 
> of the BH and location.
> In between each top of soil there are no data points in the database, 
> however, in the real world there is a soil property.  Therefore if i am 
> interpolating only the top of soil data, the interpolation method will 
> do a gradual change between the top (n) of that soil deposit and its 
> bottom which is also the top of the next underlying soil deposit (n-1).
> In order to avoid this, I would like to ensure that the 'weight' of that 
> soil deposit thickness in the interpolation is considered. Whether the 
> interpolation modules in Grass GIS can (first option)do this 
> interpolation by considering that the property extends to the next 
> underlying deposit, 


IIUC what you mean, r.to.rast3elev does exactly that: It does not do any 
interpolation. Each input map ('input' parameter) corresponds to a soil 
type. Each elevation map ('elevation' parameter) corresponds to either 
the top or the bottom level of that type. Using the -u or -l flag, you 
can ask r.to.rast3elev to fill the volume above or below your elevations 
with the value of the map.

In other words, if you have


r.mapcalc "soil1 = 1"
r.mapcalc "soil2 = 2"
r.mapcalc "soil3 = 3"

and you have 2D elevation maps resulting from the interpolation of the 
top elevation of the respective soils: elev1, elev2, elev3

You can run:

r.to.rast3elev -l input=soil1,soil2,soil3 elevation=elev1,elev2,elev3 
output=soils_3D

You should get something like this:


-----------top elevation soil1-------
1111111111111111111111111111111111111
1111111111111111111111111111111111111
1111111111111111111111111111111111111
-----------top elevation soil2-------
2222222222222222222222222222222222222
2222222222222222222222222222222222222
2222222222222222222222222222222222222
-----------top elevation soil3-------
3333333333333333333333333333333333333
3333333333333333333333333333333333333
3333333333333333333333333333333333333
3333333333333333333333333333333333333

Obviously with elevations varying in space.

There are issues if you have inversions in the order of soil layers, 
i.e. if you have something like this:

111111111111112222222222222222
111111111111112222222222222222
222222222222221111111111111111
222222222222221111111111111111

AFAIK, the only way to solve this is to cut up your area into tiles of 
homogeneous order.

Moritz