Skip to content

Commit

Permalink
adjustments and include download links
Browse files Browse the repository at this point in the history
  • Loading branch information
boriskaus committed Mar 1, 2024
1 parent 8278cb7 commit 87d31d1
Show file tree
Hide file tree
Showing 3 changed files with 45,570 additions and 47 deletions.
47 changes: 21 additions & 26 deletions docs/src/man/Tutorial_Jura.md
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ We downloaded the pdf map, and cropped it to the lower left and upper right corn
The resulting map was uploaded to zenodo; it can be downloaded with

```julia
#download_data("https://zenodo.org/10726801/SchoriM_Encl_01_Jura-map_A1.png", "SchoriM_Encl_01_Jura-map_A1.png")
download_data("https://zenodo.org/records/10726801/files/SchoriM_Encl_01_Jura-map_A1.png", "SchoriM_Encl_01_Jura-map_A1.png")
```

We also used a slightly larger version of the map along with the online tool [https://apps.automeris.io](https://apps.automeris.io) to extract the location of the corners (using the indicated blue lon/lat values on the map as reference points).
Expand Down Expand Up @@ -67,7 +67,7 @@ Within `GeophysicalModelGenerator`, we need a `longlat` coordinate system. It is
We did this and saved the resulting image in Zenodo:

```julia
#download_data("https://zenodo.org/10726801/BMes_Spline_longlat.tif", "BMes_Spline_longlat.tif")
download_data("https://zenodo.org/records/10726801/files/BMes_Spline_longlat.tif", "BMes_Spline_longlat.tif")
```

Now, import the GeoTIFF as:
Expand All @@ -86,7 +86,7 @@ Often that is not the case and you have to use the mapview along with the digiti
As example, we use the cross-section

```julia
#download_data("https://zenodo.org/10726801/BMes_Spline_longlat.tif", "BMes_Spline_longlat.tif")
download_data("https://zenodo.org/records/10726801/files/Schori_2020_Ornans-Miserey-v2_whiteBG.png", "Schori_2020_Ornans-Miserey-v2_whiteBG.png")
Corner_LowerLeft = (5.92507, 47.31300, -2.0)
Corner_UpperRight = (6.25845, 46.99550, 2.0)
CrossSection_1 = Screenshot_To_GeoData("Schori_2020_Ornans-Miserey-v2_whiteBG.png", Corner_LowerLeft, Corner_UpperRight) # name should have "colors" in it
Expand Down Expand Up @@ -202,15 +202,15 @@ Write_Paraview(CrossSection_1_cart,"CrossSection_1_cart")
The result looks like:
![Jura_Tutorial_1](../assets/img/Jura_1.png)

## 3. Geological block model
Yet, if you want to perform a numerical simulation of the Jura, it is more convenient to rotate the maps such that we can perform a simulation perpendicular to the strike of the mountain belt.
This can be done with `RotateTranslateScale`:

```julia
RotationAngle = -43;
TopoGeology_cart_rot = RotateTranslateScale(TopoGeology_cart, Rotate=RotationAngle);
Basement_cart_rot = RotateTranslateScale(Basement_cart, Rotate=RotationAngle);
CrossSection_1_cart_rot = RotateTranslateScale(CrossSection_1_cart, Rotate=RotationAngle);
nothing #hide
RotationAngle = -43
TopoGeology_cart_rot = RotateTranslateScale(TopoGeology_cart, Rotate=RotationAngle)
Basement_cart_rot = RotateTranslateScale(Basement_cart, Rotate=RotationAngle)
CrossSection_1_cart_rot = RotateTranslateScale(CrossSection_1_cart, Rotate=RotationAngle)
```

Next, we can create a new computational grid that is more conveniently oriented:
Expand All @@ -219,9 +219,8 @@ We create both a surface and a 3D block
```julia
nx, ny, nz = 1024, 1024, 128
x,y,z = range(-100,180,nx), range(-50,70,ny), range(-8,4,nz)
ComputationalSurf = CartData(XYZGrid(x,y,0));
ComputationalGrid = CartData(XYZGrid(x,y,z));
nothing #hide
ComputationalSurf = CartData(XYZGrid(x,y,0))
ComputationalGrid = CartData(XYZGrid(x,y,z))
```

Re-interpolate the rotated to the new grid:
Expand All @@ -235,24 +234,21 @@ Next we can use the surfaces to create a 3D block model.
We start with a block model that has the different rocktypes:

```julia
Phases = zeros(Int8,size(ComputationalGrid.x)); #Define rock types
nothing #hide
Phases = zeros(Int8,size(ComputationalGrid.x)) #Define rock types
```

Set everything below the topography to 1

```julia
id = BelowSurface(ComputationalGrid, GeologyTopo_comp_surf);
Phases[id] .= 1;
nothing #hide
id = BelowSurface(ComputationalGrid, GeologyTopo_comp_surf)
Phases[id] .= 1
```

The basement is set to 2

```julia
id = BelowSurface(ComputationalGrid, Basement_comp_surf);
Phases[id] .= 2;
nothing #hide
id = BelowSurface(ComputationalGrid, Basement_comp_surf)
Phases[id] .= 2
```

Add to the computational grid:
Expand All @@ -265,18 +261,17 @@ ComputationalGrid = RemoveField(ComputationalGrid,"Z")
Save the surfaces, cross-section and the grid:

```julia
Write_Paraview(GeologyTopo_comp_surf,"GeologyTopo_comp_surf");
Write_Paraview(Basement_comp_surf, "Basement_comp_surf");
Write_Paraview(CrossSection_1_cart_rot,"CrossSection_1_cart_rot");
Write_Paraview(ComputationalGrid,"ComputationalGrid");
nothing #hide
Write_Paraview(GeologyTopo_comp_surf,"GeologyTopo_comp_surf")
Write_Paraview(Basement_comp_surf, "Basement_comp_surf")
Write_Paraview(CrossSection_1_cart_rot,"CrossSection_1_cart_rot")
Write_Paraview(ComputationalGrid,"ComputationalGrid")
```

We can visualize this in paraview:
![Jura_Tutorial_2](../assets/img/Jura_2.png)

Note that the `y`-direction is now perpendicular to the Jura mountains.
The paraview statefiles to generate the two figures shown here are `/tutorials/Jura_1.pvsm` and `/tutorials/Jura_2.pvsm`.
We use a vertical exaggeration of factor two. Also note that the `y`-direction is now perpendicular to the Jura mountains.
The paraview statefiles to generate this figure is `/tutorials/Jura_2.pvsm`.

---

Expand Down
Loading

0 comments on commit 87d31d1

Please sign in to comment.