Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
/test/*.vts
/test/*.vtu
/.vscode
.DS_Store
Manifest.toml
/docs/build
Expand Down
6 changes: 3 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -50,18 +50,18 @@ Downloads = "1"
FFMPEG = "0.4"
FileIO = "1"
GDAL_jll = "300.900.0 - 301.901.0"
GLMakie = "0.8, 0.9, 0.10, 0.11"
GLMakie = "0.8 - 0.11"
GeoParams = "0.2 - 0.7"
Geodesy = "1"
GeometryBasics = "0.1 - 0.4"
GeometryBasics = "0.1 - 0.5"
Glob = "1.2 - 1.3"
GMT = "1"
GridapGmsh = "0.5 - 0.7"
ImageIO = "0.1 - 0.6"
Interpolations = "0.14, 0.15"
JLD2 = "0.4, 0.5"
LightXML = "0.8, 0.9"
MeshIO = "0.1 - 0.4"
MeshIO = "0.1 - 0.5"
NCDatasets = "0.14"
NearestNeighbors = "0.2 - 0.4"
Parameters = "0.9 - 0.12"
Expand Down
12 changes: 6 additions & 6 deletions docs/src/man/Tutorial_Basic.md
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ GeoData
````

This is a so-called `GeoData` object, which is a 3D grid of seismic velocities as a function of longitude, latitude and depth, which can include various fields (here we only have a single field: `:dVp_Percentage`)
We can save this in `VTK` format, which is a widely used format that can for exampke be read by the 3D open-source visualization tool [Paraview](https://www.paraview.org/):
We can save this in `VTK` format, which is a widely used format that can for example be read by the 3D open-source visualization tool [Paraview](https://www.paraview.org/):

```julia
write_paraview(Tomo_Alps_full,"Tomo_Alps_full")
Expand Down Expand Up @@ -101,7 +101,7 @@ After loading the new data again in paraview, switching to the proper data field
![Basic_Tutorial_2](../assets/img/Basic_Tutorial_2.png)

### 3. Create cross sections
Paraview has the option to `Slice` through the data but it is not very intuitive to do this in 3D. Another limitation of Paraview is that it does not have native support for spherical coordinates, and therefore the data is translated to cartesian (`x`,`y`,`z`) coordinates (with the center of the Earth at `(0,0,0)`).
Paraview has the option to `Slice` through the data but it is not very intuitive to do this in 3D. Another limitation of Paraview is that it does not have native support for spherical coordinates, and therefore the data is translated to Cartesian (`x`,`y`,`z`) coordinates (with the center of the Earth at `(0,0,0)`).
That makes this a bit cumbersome to make a cross-section at a particular location.
If you are interested in this you can use the `cross_section` function:

Expand Down Expand Up @@ -214,9 +214,9 @@ After doing all these steps, you should see something like this:
![Basic_Tutorial_3](../assets/img/Basic_Tutorial_3.png)

### 4. Cartesian data
As you can see, the curvature or the Earth is taken into account here. Yet, for many applications it is more convenient to work in Cartesian coordinates (kilometers) rather then in geographic coordinates.
As you can see, the curvature or the Earth is taken into account here. Yet, for many applications it is more convenient to work in Cartesian coordinates (kilometers) rather than in geographic coordinates.
`GeophysicalModelGenerator` has a number of tools for this.
First we need do define a `ProjectionPoint` around which we project the data
First we need to define a `ProjectionPoint` around which we project the data

```julia
proj = ProjectionPoint(Lon=12.0,Lat =43)
Expand Down Expand Up @@ -267,7 +267,7 @@ Saved file: Topo_cart.vts
As the coordinates are now aligned with the `x`,`y`,`z` coordinate axes in Paraview it is now straightforward to use the build-in tools to explore the data.

### 5. Rectilinear data
Yet, because of the curvature of the Earth, the resulting 3D model is not strictly rectilinear, which is often a requiment for cartesian numerical models.
Yet, because of the curvature of the Earth, the resulting 3D model is not strictly rectilinear, which is often a requirement for Cartesian numerical models.
This can be achieved in a relatively straightforward manner, by creating a new 3D dataset that is slightly within the curved boundaries of the projected data set:

```julia
Expand Down Expand Up @@ -322,7 +322,7 @@ Saved file: Topo_rect.vts

![Basic_Tutorial_5](../assets/img/Basic_Tutorial_5.png)

At this stage, the data can directly be used to generate cartesian numerical model setups, as explained in the other tutorials.
At this stage, the data can directly be used to generate Cartesian numerical model setups, as explained in the other tutorials.

---

Expand Down
6 changes: 3 additions & 3 deletions docs/src/man/Tutorial_Jura.md
Original file line number Diff line number Diff line change
Expand Up @@ -94,8 +94,8 @@ CrossSection_1 = screenshot_to_GeoData("Schori_2020_Ornans-Miserey-v2_whiteBG.pn

Note that we slightly modified the image to save it with a white instead of a transparent background

# 2. Project the data to a cartesian grid
At this stage, we have all data in geographic coordinates. In most cases it is more useful to have them in cartesian coordinates.
# 2. Project the data to a Cartesian grid
At this stage, we have all data in geographic coordinates. In most cases it is more useful to have them in Cartesian coordinates.
Moreover, the resolution of the grids is different. Whereas the `TopoGeology` has a size of `(3721, 2641, 1)`, `Basement` has size `(2020, 1751, 1)`.
It is often useful to have them on exactly the same size grid

Expand Down Expand Up @@ -169,7 +169,7 @@ CartData
fields : (:Basement,)
```

Finally, we can also transfer the cross-section to cartesian coordinates. As this is just for visualization, we will
Finally, we can also transfer the cross-section to Cartesian coordinates. As this is just for visualization, we will
use `convert2CartData` in this case

```julia
Expand Down
6 changes: 3 additions & 3 deletions docs/src/man/Tutorial_LaPalma.md
Original file line number Diff line number Diff line change
Expand Up @@ -60,22 +60,22 @@ write_paraview(Topo,"Topo")
![LaPalma_EQTopo_GeoData](../assets/img/TopoEQs_LaPalma_GeoData.png)
Note that this data is in geographic coordinates, which makes it non-trivial to create slices through the data (see coordinate axis in the plot, where `z` is *not* pointing upwards).

## 2. Convert data to cartesian coordinates
## 2. Convert data to Cartesian coordinates
In order to create model setups, it is helpful to first transfer the data to Cartesian.
This requires us to first determine a *projection point*, that is fixed. Often, it is helpful to use the center of the topography for this. In the present example, we will center the model around La Palma itself:

```julia
proj = ProjectionPoint(Lon=-17.84, Lat=28.56)
```

Once this is done you can convert the topographic data to the cartesian reference frame
Once this is done you can convert the topographic data to the Cartesian reference frame

```julia
EQ_cart = convert2CartData(data_all_EQ, proj);
Topo_cart = convert2CartData(Topo, proj)
```

It is important to realize that the cartesian coordinates of the topographic grid is no longer strictly orthogonal after this conversion. You don't notice that in the current example, as the model domain is rather small.
It is important to realize that the Cartesian coordinates of the topographic grid is no longer strictly orthogonal after this conversion. You don't notice that in the current example, as the model domain is rather small.
In other cases, however, this is quite substantial (e.g., India-Asia collision zone).
LaMEM needs an orthogonal grid of topography, which we can create with:

Expand Down
4 changes: 2 additions & 2 deletions docs/src/man/Tutorial_NumericalModel_2D.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ The aim of this tutorial is to show you how to create 2D numerical model setups

### 2D Subduction setup

Lets start with creating a 2D model setup in cartesian coordinates, which uses the `CartData` data structure
Lets start with creating a 2D model setup in Cartesian coordinates, which uses the `CartData` data structure

```julia
using GeophysicalModelGenerator
Expand Down Expand Up @@ -351,7 +351,7 @@ We have a number of other functions to help create a geometry, specifically:
- `AddCylinder!`

The help functions are quite self-explanatory, so we won't show it in detail here.
If you have a topography surface or any other horizontal surface, you can surface with the cartesian grid with `above_surface` or `below_surface`.
If you have a topography surface or any other horizontal surface, you can surface with the Cartesian grid with `above_surface` or `below_surface`.

Also, if you wish to take a seismic tomography as inspiration to set a slab geometry, you can interpolate it to a `CartGrid` with the same dimensions and use that with the julia `findall` function.

Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/Tutorial_NumericalModel_3D.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ The aim of this tutorial is to show you how to create 3D numerical model setups

### 3D Subduction setup

Lets start with creating a 3D model setup in cartesian coordinates, which uses the `CartData` data structure
Lets start with creating a 3D model setup in Cartesian coordinates, which uses the `CartData` data structure

```julia
using GeophysicalModelGenerator
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/Tutorial_VolcanoModel_3D.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ The aim of this tutorial is to show you how to create 3D numerical model setups

### Generating the model

Lets start with creating a 3D model setup in cartesian coordinates, which uses the `CartData` data structure, with a resolution of $ 128 \times 128 \times 128 $ grid points, inside the domain $\Omega \in [-100,100] \times [-100,100] \times [-110,50]$ km
Lets start with creating a 3D model setup in Cartesian coordinates, which uses the `CartData` data structure, with a resolution of $ 128 \times 128 \times 128 $ grid points, inside the domain $\Omega \in [-100,100] \times [-100,100] \times [-110,50]$ km

```julia
using GeophysicalModelGenerator
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/datastructures.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

The main data structure used in GeophysicalModelGenerator.jl is `GeoData`, which contains info about the `longitude`,`latitude`, and `depth` of a data set, as well as several data sets itself.
We also provide a `UTMData`, which is essentially the same but with UTM coordinates, and a `CartData` structure, which has Cartesian coordinates in kilometers (as used in many geodynamic codes). If one wishes to transfer `GeoData` to `CartData`, one needs to provide a `ProjectionPoint`.
For plotting, we transfer this into the `ParaviewData` structure, which has cartesian coordinates centered around the center of the Earth. We employ the `wgs84` reference ellipsoid as provided by the [Geodesy.jl](https://github.com/JuliaGeo/Geodesy.jl) package to perform this transformation.
For plotting, we transfer this into the `ParaviewData` structure, which has Cartesian coordinates centered around the center of the Earth. We employ the `wgs84` reference ellipsoid as provided by the [Geodesy.jl](https://github.com/JuliaGeo/Geodesy.jl) package to perform this transformation.

```@docs
GeoData
Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/lamem.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# LaMEM

In order to generate geodynamic simulations from setups created with `GeophysicalModelGenerator.jl`, we provide a few routines that directly create marker input files for the 3D geodynamic modelling software [LaMEM](https://github.com/UniMainzGeo/LaMEM), which is an open-source cartesian code to perform crustal and lithospheric-scale simulations.
In order to generate geodynamic simulations from setups created with `GeophysicalModelGenerator.jl`, we provide a few routines that directly create marker input files for the 3D geodynamic modelling software [LaMEM](https://github.com/UniMainzGeo/LaMEM), which is an open-source Cartesian code to perform crustal and lithospheric-scale simulations.
If you want to learn how to run LaMEM simulations, the easiest way to get started is by looking at [LaMEM.jl](https://github.com/JuliaGeodynamics/LaMEM.jl) which is integrated with `GMG`

The routines provided here have the following functionality:
Expand Down
6 changes: 3 additions & 3 deletions docs/src/man/projection.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ Typically, you load a dataset by reading it into julia and either generating a `

If you write the data to `Paraview`, it is internally converted to a Paraview structure (which involves `x,y,z` Cartesian Earth-Centered-Earth-Fixed (ECEF) coordinates using the `wgs84` ellipsoid).

Yet, if you do geodynamic calculations the chances are that the geodynamic code does not operate in spherical coordinates, but rather use cartesian ones. In that case you should transfer your data to the `CartData` structure, which requires you to specify a `ProjectionPoint` that is a point on the map that will later have the coordinates `(0,0)` in the `CartData` structure.
Yet, if you do geodynamic calculations the chances are that the geodynamic code does not operate in spherical coordinates, but rather use Cartesian ones. In that case you should transfer your data to the `CartData` structure, which requires you to specify a `ProjectionPoint` that is a point on the map that will later have the coordinates `(0,0)` in the `CartData` structure.


#### 1. Converting
Expand Down Expand Up @@ -75,7 +75,7 @@ This shows that the model is ~5600 by 3000 km.
Whereas this is ok to look at and compare with a LaMEM model setup, we cannot use it to perform internal calculations (or to generate a LaMEM model setup), because the `x` and `y` coordinates are distorted and not orthogonal.

#### 2. Projecting data
For use with LaMEM, you would need an orthogonal cartesian grid. From the last command above we get some idea on the area, so we can create this:
For use with LaMEM, you would need an orthogonal Cartesian grid. From the last command above we get some idea on the area, so we can create this:
```julia-repl
julia> Topo_Cart_orth = CartData(xyz_grid(-2000:20:2000,-1000:20:1000,0))
CartData
Expand All @@ -97,7 +97,7 @@ CartData
julia> write_paraview(Topo_Cart_orth,"Topo_Cart_orth");
```
![Topo_Europe_CartData_Proj](../assets/img/Topo_Europe_CartData_Proj.png)
So this interpolates the topographic data from the `GeoData` to the orthogonal cartesian grid (which can be used with LaMEM, for example).
So this interpolates the topographic data from the `GeoData` to the orthogonal Cartesian grid (which can be used with LaMEM, for example).

You can do similar projections with full 3D data sets or pointwise data.

Expand Down
2 changes: 1 addition & 1 deletion docs/src/man/tutorial_load3DSeismicData.md
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,7 @@ If you want to plot iso-surfaces (e.g. at 3%), you can use the `Clip` option aga

![Paraview_6](../assets/img/Tutorial_Zhao_Paraview_6.png)

Note that using geographic coordinates is slightly cumbersome in Paraview. If your area is not too large, it may be advantageous to transfer all data to cartesian coordinates, in which case it is easier to create slices through the model. This is explained in some of the other tutorials.
Note that using geographic coordinates is slightly cumbersome in Paraview. If your area is not too large, it may be advantageous to transfer all data to Cartesian coordinates, in which case it is easier to create slices through the model. This is explained in some of the other tutorials.

#### 5. Extract and plot cross-sections of the data
In many cases you would like to create cross-sections through the 3D data sets as well, and visualize that in Paraview. That is in principle possible in Paraview as well (using the `Slice` tool, as described above). Yet, in many cases we want to have it at a specific depth, or through pre-defined `lon/lat` coordinates.
Expand Down
12 changes: 6 additions & 6 deletions docs/src/man/tutorial_local_Flegrei.md
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# Campi Flegrei Volcano tutorial using cartesian coordinates
# Campi Flegrei Volcano tutorial using Cartesian coordinates

## Goal

This tutorial visualizes available 3D data at a local volcano (Campi Flegrei caldera, Italy) using cartesian coordinates. This is done, e.g., when we only have geomorphological data in cartesian coordinates. It includes geological and geophysical data in UTM format from the following papers:
This tutorial visualizes available 3D data at a local volcano (Campi Flegrei caldera, Italy) using Cartesian coordinates. This is done, e.g., when we only have geomorphological data in Cartesian coordinates. It includes geological and geophysical data in UTM format from the following papers:

- Two shape files containing coastline and faults:
- Vilardo, G., Ventura, G., Bellucci Sessa, E. and Terranova, C., 2013. Morphometry of the Campi Flegrei caldera (southern Italy). Journal of maps, 9(4), pp.635-640. doi:10.1080/17445647.2013.842508
Expand All @@ -29,11 +29,11 @@ Make sure that you are in the unzipped directory.

#### 2. Geomorphology

Load both the the shape (.shp) files contained in `./Geomorphology/*.shp` inside Paraview. In the following figures we show the Cartesian representation (not geolocalized - left) and the UTM (UTM). Our shape files can only be loaded in the cartesian:
Load both the the shape (.shp) files contained in `./Geomorphology/*.shp` inside Paraview. In the following figures we show the Cartesian representation (not geolocalized - left) and the UTM (UTM). Our shape files can only be loaded in the Cartesian:

![Tutorial_Flegrei_Geomorphology](../assets/img/Flegrei_Geomorphology.png)

To reproduce it, represent the coastline as data points with black solid color and assign your favourite color map to the morphology. Note that each block color number corresponds to a different morphology. Beware that this file only works in cartesian coordinate, as it is still impossible to generate shape files in real UTM coordinates
To reproduce it, represent the coastline as data points with black solid color and assign your favourite color map to the morphology. Note that each block color number corresponds to a different morphology. Beware that this file only works in Cartesian coordinate, as it is still impossible to generate shape files in real UTM coordinates

#### 3. Earthquakes

Expand All @@ -55,15 +55,15 @@ julia> EQ_Data_UTM = UTMData(WE, SN, depth, 33, true, (Depth=depth * m,T
julia> Data_set_UTM = convert(GeophysicalModelGenerator.GeoData,EQ_Data_UTM)
julia> write_paraview(Data_set_UTM, "CF_Earthquakes_UTM", PointsData=true)
```
Save in paraview with both cartesian and UTM formats. The final seismicity map looks like this:
Save in paraview with both Cartesian and UTM formats. The final seismicity map looks like this:

![Tutorial_Flegrei_seismicity](../assets/img/Flegrei_Seismicity.png)

The colour scale distinguishes earthquakes of different decades. Notice the progressive migration of recent seismicity (black dots) towards East.

#### 4. Velocity model

Using the Alps tutorial it is easy to create a paraview file from the Vp, Vs and Vp/Vs model in `./TravelTmeTomography/modvPS.dat` for both cartesian and UTM coordinates.
Using the Alps tutorial it is easy to create a paraview file from the Vp, Vs and Vp/Vs model in `./TravelTmeTomography/modvPS.dat` for both Cartesian and UTM coordinates.

```julia
julia> using DelimitedFiles, GeophysicalModelGenerator
Expand Down
4 changes: 2 additions & 2 deletions docs/src/man/visualise.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ julia> using GeophysicalModelGenerator, GLMakie
julia> using GMT, JLD2
```
The last line loads packages we need to read in the pre-generated `JLD2` file (see the tutorials) and download topography from the region.
Lets load the data, by first moving to the correct directory (will likely be different on your machine).
Let's load the data, by first moving to the correct directory (will likely be different on your machine).
```julia
julia> ;
shell> cd ~/Downloads/Zhao_etal_2016_data/
Expand All @@ -34,7 +34,7 @@ Which will look like:
![Tutorial_Visualize](../assets/img/Tutorial_Visualize.png)


This is an example where we used `GeoData` to visualize results. Alternatively, we can also visualize results in km (often more useful for numerical modelling setups). For `Visualize` to work with this, we however need orthogonal cartesian data, which can be obtained by projecting both the data.
This is an example where we used `GeoData` to visualize results. Alternatively, we can also visualize results in km (often more useful for numerical modelling setups). For `Visualize` to work with this, we however need orthogonal Cartesian data, which can be obtained by projecting both the data.
```julia
julia> p=ProjectionPoint(Lon=10, Lat=45)
ProjectionPoint(45.0, 10.0, 578815.302916711, 4.983436768349297e6, 32, true)
Expand Down
5 changes: 4 additions & 1 deletion src/stl.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,10 @@
# These are routines that allow importing *.stl triangular surfaces and employ them

using FileIO
using GeometryBasics: TriangleP, Mesh, normals, PointMeta, coordinates
# using GeometryBasics: TriangleP, Mesh, normals, PointMeta, coordinates
# I removed TriangleP and PointMeta, as they are not existing in GeometryBasics anymore (past 0.5.0)
# HD, 2024-06-10
using GeometryBasics: Mesh, normals, coordinates
using LinearAlgebra

# Warning: the TriangleIntersect dependency does not seem to work on different machines, as the developer did not add a version number..
Expand Down
Loading