diff --git a/.gitignore b/.gitignore index 157136df..91c96fe2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ /test/*.vts /test/*.vtu +/.vscode .DS_Store Manifest.toml /docs/build diff --git a/Project.toml b/Project.toml index 3ef0224d..8361ec0e 100644 --- a/Project.toml +++ b/Project.toml @@ -50,10 +50,10 @@ 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" @@ -61,7 +61,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" diff --git a/docs/src/man/Tutorial_Basic.md b/docs/src/man/Tutorial_Basic.md index 958e55c0..52f30cb0 100644 --- a/docs/src/man/Tutorial_Basic.md +++ b/docs/src/man/Tutorial_Basic.md @@ -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") @@ -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: @@ -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) @@ -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 @@ -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. --- diff --git a/docs/src/man/Tutorial_Jura.md b/docs/src/man/Tutorial_Jura.md index 8e629457..dbc49130 100644 --- a/docs/src/man/Tutorial_Jura.md +++ b/docs/src/man/Tutorial_Jura.md @@ -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 @@ -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 diff --git a/docs/src/man/Tutorial_LaPalma.md b/docs/src/man/Tutorial_LaPalma.md index 5ca85168..018c8e1e 100644 --- a/docs/src/man/Tutorial_LaPalma.md +++ b/docs/src/man/Tutorial_LaPalma.md @@ -60,7 +60,7 @@ 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: @@ -68,14 +68,14 @@ This requires us to first determine a *projection point*, that is fixed. Often, 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: diff --git a/docs/src/man/Tutorial_NumericalModel_2D.md b/docs/src/man/Tutorial_NumericalModel_2D.md index 6e767171..cd7b560c 100644 --- a/docs/src/man/Tutorial_NumericalModel_2D.md +++ b/docs/src/man/Tutorial_NumericalModel_2D.md @@ -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 @@ -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. diff --git a/docs/src/man/Tutorial_NumericalModel_3D.md b/docs/src/man/Tutorial_NumericalModel_3D.md index 4e26e292..6820b2d3 100644 --- a/docs/src/man/Tutorial_NumericalModel_3D.md +++ b/docs/src/man/Tutorial_NumericalModel_3D.md @@ -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 diff --git a/docs/src/man/Tutorial_VolcanoModel_3D.md b/docs/src/man/Tutorial_VolcanoModel_3D.md index aaf05a2f..3fd594ce 100644 --- a/docs/src/man/Tutorial_VolcanoModel_3D.md +++ b/docs/src/man/Tutorial_VolcanoModel_3D.md @@ -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 diff --git a/docs/src/man/datastructures.md b/docs/src/man/datastructures.md index 5b170ea1..eb7576e0 100644 --- a/docs/src/man/datastructures.md +++ b/docs/src/man/datastructures.md @@ -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 diff --git a/docs/src/man/lamem.md b/docs/src/man/lamem.md index 7f081051..ac5a909c 100644 --- a/docs/src/man/lamem.md +++ b/docs/src/man/lamem.md @@ -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: diff --git a/docs/src/man/projection.md b/docs/src/man/projection.md index f8934ffd..4cb74de4 100644 --- a/docs/src/man/projection.md +++ b/docs/src/man/projection.md @@ -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 @@ -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 @@ -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. diff --git a/docs/src/man/tutorial_load3DSeismicData.md b/docs/src/man/tutorial_load3DSeismicData.md index 12a5b0d4..a40cbe54 100644 --- a/docs/src/man/tutorial_load3DSeismicData.md +++ b/docs/src/man/tutorial_load3DSeismicData.md @@ -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. diff --git a/docs/src/man/tutorial_local_Flegrei.md b/docs/src/man/tutorial_local_Flegrei.md index fb3adb68..655a2881 100644 --- a/docs/src/man/tutorial_local_Flegrei.md +++ b/docs/src/man/tutorial_local_Flegrei.md @@ -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 @@ -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 @@ -55,7 +55,7 @@ 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) @@ -63,7 +63,7 @@ The colour scale distinguishes earthquakes of different decades. Notice the prog #### 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 diff --git a/docs/src/man/visualise.md b/docs/src/man/visualise.md index 4f5cd9ac..e9524256 100644 --- a/docs/src/man/visualise.md +++ b/docs/src/man/visualise.md @@ -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/ @@ -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) diff --git a/src/stl.jl b/src/stl.jl index 6762e60d..2a24bafb 100644 --- a/src/stl.jl +++ b/src/stl.jl @@ -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..