Hi Moritz, Thanks for your email. What you described is the concept for developing most geological model by correlating top of layers (rock, soils, formation etc.) which is the correct approach when geological contacts are clear.
In my case (as part of a masters/engineering degree), I have so much data and i am interpolating soils, that i am trying another approach where i will be interpolating a property in 3D. Although there are some regional conceptual models in my areas, we often find that reality differs quite a bit from the conceptual understanding of the area. Therefore i am attempting to develop a ''newish" approach by interpolating the data set by using a numerical property for the soil type.
I will be interpolating most likely particle size (x mm for sand, x mm for silts, x mm for clays etc), and hydraulic conductivity (K). This could be any parameters such as moisture content, geotech blow counts, concentrations etc.
That said, from the borehole information, the soil type or particle size is constant until the next underlying formation. However, the data base only provide the property at the top of the formation (ex: silt) like typical BH databases. Once the interpolation is completed, i can identify the soils by the color coding and custom scales.
So i am wondering if in Grass, i need to generate a property (particle size, K) for each meter of elevation so that i have a data point in between the top and bottom of the layer, or if there are some modules that can automatically consider that the property extends to the next underlying layer.
I think from what i understand, i will need to generate a data point for each meter of elevation along the borehole axis, and possibly interpolate each elevation separately or using a module that can do 3D interpolation such as v.vol.rst. although i want to constrain the interpolation horizontally, so I may need to interpolate each elevation separately (manually or in batch mode).
alternatively, i could be able to do a query in grass at each elevation to search for the property above along the Bh axis during the interpolation process but this seems complicated. I can generate properties for each elevations in Access in import into Grass.
Thanks
Francois
Le mar. 10 juil. 2018 à 21:59, Moritz Lennert <mlennert@club.worldonline.be> a écrit :
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
