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
_______________________________________________
grass-user mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/grass-user
