Add surfaces routines

AUTHORS.md

@@ -1,6 +1,6 @@
 # Authors
 
-GeophysicalModelGenerator.jl's development is coordinated by a group of *principal developers*,
+`GeophysicalModelGenerator.jl`'s development is coordinated by a group of *principal developers*,
 who are also its main contributors and who can be contacted in case of
 questions about GeophysicalModelGenerator.jl. In addition, there are *contributors* who have
 provided substantial additions or modifications. Together, these two groups form
@@ -14,11 +14,14 @@ provided substantial additions or modifications. Together, these two groups form
 
 
 ## Contributors
-The following people contributed major additions or modifications to GeophysicalModelGenerator.jl and
+The following people contributed major additions or modifications to `GeophysicalModelGenerator.jl` and
 are listed in alphabetical order:
 
 * Pascal Aellig
+* Albert De Montserrat
 * Luca De Siena
+* Jacob Frasunkiewicz
+* Lukas Fuchs
 * Andrea Piccolo
 * Hendrik Rachona
 * Christian Schuler

docs/make.jl

@@ -84,7 +84,7 @@ makedocs(;
             "Overview" =>  "man/tutorials.md",
             "1 - 3D seismic tomography from ASCII" =>  "man/tutorial_load3DSeismicData.md",
             "2 - 3D seismic tomography from netCDF" =>  "man/tutorial_loadregular3DSeismicData_netCDF.md",
-            "3 - Visualize Moho topography" =>  "man/tutorial_MohoTopo.md",
+            "3 - Visualize Moho topography" =>  "man/Tutorial_MohoTopo_Spada.md",
             "4 - Create GMT-based topography" =>  "man/tutorial_GMT_Topography.md",
             "5 - Coastlines" =>  "man/tutorial_Coastlines.md",
             "6 - Import screenshots" =>  "man/tutorial_Screenshot_To_Paraview.md",
@@ -97,14 +97,15 @@ makedocs(;
             "13 - Campi Flegrei" =>  "man/tutorial_local_Flegrei.md",
             "14 - La Palma volcano Model" =>  "man/Tutorial_LaPalma.md",
             "15 - Create movies" =>  "man/tutorial_time_Seismicity.md",
-            "16 - Fault Density Map" =>  "man/tutorial_Fault_Map.md",
+            "16 - Fault Density Map" =>  "man/Tutorial_FaultDensity.md",
             "17 - Jura tutorial" =>  "man/Tutorial_Jura.md"
         ],
         "User Guide" => Any[
             "Installation" =>  "man/installation.md",
             "Data Structures" =>  "man/datastructures.md",
             "Data Import" =>  "man/dataimport.md",
             "Projection" =>  "man/projection.md",
+            "Surfaces" =>  "man/surfaces.md",
             "Paraview output" => "man/paraview_output.md",
             "Tools" => "man/tools.md",
             "Visualisation" => "man/visualise.md",

docs/src/man/Tutorial_FaultDensity.md

@@ -0,0 +1,114 @@
+```@meta
+EditURL = "../../../tutorials/Tutorial_FaultDensity.jl"
+```
+
+# Fault Density Map
+
+## Goal
+In this tutorial, Fault Data is loaded as Shapefiles, which is then transformed to raster data. With the help of that a fault density map of Europe is created.
+
+## 1. Load Data
+
+Load packages:
+
+```julia
+using GeophysicalModelGenerator, Shapefile, Plots, Rasters, GeoDatasets, Interpolations
+```
+
+Data is taken from "Active Faults of Eurasia Database AFEAD v2022" DOI:10.13140/RG.2.2.25509.58084
+You need to download it manually (as it doesn't seem to work automatically), and can load the following file:
+
+```julia
+File   = "AFEAD_v2022/AFEAD_v2022.shp"
+```
+
+Load data using the `Shapefile` package:
+
+```julia
+table  = Shapefile.Table(File)
+geoms  = Shapefile.shapes(table)
+CONF   = table.CONF
+```
+
+Raster the shapefile data
+
+```julia
+ind    = findall((table.CONF .== "A") .| (table.CONF .== "B") .| (table.CONF .== "C"))
+faults = Shapefile.Handle(File).shapes[ind]
+faults = rasterize(last,faults; res=(0.12,0.12), missingval=0, fill=1, atol = 0.4, shape=:line)
+lon    = faults.dims[1]
+lat    = faults.dims[2]
+```
+
+Download coastlines with `GeoDatasets`:
+
+```julia
+lonC,latC,dataC = GeoDatasets.landseamask(;resolution='l',grid=10)
+```
+
+Interpolate to fault grid
+
+```julia
+itp        = linear_interpolation((lonC, latC), dataC)
+coastlines = itp[lon.val,lat.val]
+coastlines = map(y -> y > 1 ? 1 : y, coastlines)
+```
+
+Plot the fault data
+
+```julia
+heatmap(lon.val,lat.val,coastlines',legend=false,colormap=cgrad(:gray1,rev=true),alpha=0.4);
+plot!(faults; color=:red,legend = false,title="Fault Map World",ylabel="Lat",xlabel="Lon")
+```
+
+![tutorial_Fault_Map](../assets/img/WorldMap.png)
+
+Restrict area to Europe
+
+```julia
+indlat = findall((lat .> 35) .& (lat .< 60))
+Lat    = lat[indlat]
+indlon = findall((lon .> -10) .& (lon .< 35))
+Lon    = lon[indlon]
+data   = faults.data[indlon,indlat]
+```
+
+Create GeoData from restricted data
+
+```julia
+Lon3D,Lat3D, Faults = LonLatDepthGrid(Lon,Lat,0);
+Faults[:,:,1]       = data
+Data_Faults         = GeoData(Lon3D,Lat3D,Faults,(Faults=Faults,))
+```
+
+## 2. Create Density Map
+Create a density map of the fault data. This is done with the `CountMap` function. This function takes a specified field of a 2D `GeoData` struct and counts the entries in all control areas which are defined by steplon (number of control areas in lon direction) and steplat (number of control areas in lat direction). The field should only consist of 0.0 and 1.0 and the steplength. The final result is normalized by the highest count.
+
+```julia
+steplon  = 188
+steplat  = 104
+cntmap   = CountMap(Data_Faults,"Faults",steplon,steplat)
+```
+
+Plot the density map with coastlines
+
+```julia
+lon = unique(cntmap.lon.val)
+lat = unique(cntmap.lat.val)
+coastlinesEurope = itp[lon,lat]
+coastlinesEurope = map(y -> y > 1 ? 1 : y, coastlinesEurope)
+```
+
+Plot this using `Plots.jl`:
+
+```julia
+heatmap(lon,lat,coastlinesEurope',colormap=cgrad(:gray1,rev=true),alpha=1.0);
+heatmap!(lon,lat,cntmap.fields.CountMap[:,:,1]',colormap=cgrad(:batlowW,rev=true),alpha = 0.8,legend=true,title="Fault Density Map Europe",ylabel="Lat",xlabel="Lon")
+```
+
+![tutorial_Fault_Map](../assets/img/FaultDensity.png)
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/man/Tutorial_Jura.md

@@ -2,7 +2,7 @@
 EditURL = "../../../tutorials/Tutorial_Jura.jl"
 ```
 
-# 17 - Create a model for the Jura mountains
+# Create a 3D model of the Jura mountains
 
 ## Aim
 In this tutorial, your will learn how to use drape a geological map on top of a digital topography model, import GeoTIFF surfaces and add cross-sections from screenshots to the model setup.
@@ -47,7 +47,7 @@ Geology  = Screenshot_To_GeoData("SchoriM_Encl_01_Jura-map_A1.png", lowerleft, u
 You can "drape" this image on the topographic map with
 
 ```julia
-TopoGeology = DrapeOnTopo(Topo, Geology)
+TopoGeology = drape_on_topo(Topo, Geology)
 ```
 
 ```julia
@@ -77,7 +77,7 @@ Basement = ImportGeoTIFF("BMes_Spline_longlat.tif", fieldname=:Basement, removeN
 ```
 
 the `removeNaN_z` option removes `NaN` values from the dataset and instead uses the z-value of the nearest point.
-That is important if you want to use this surface to generate a 3D model setup (using `BelowSurface`, for example).
+That is important if you want to use this surface to generate a 3D model setup (using `belowSurface`, for example).
 
 The thesis also provides a few interpreted vertical cross-sections. As before, we import them as a screenshot and estimate the lower-left and upper right corners.
 In this particular case, we are lucky that the `lon/lat` values are indicated on the cross-section.
@@ -177,10 +177,10 @@ CrossSection_1_cart = Convert2CartData(CrossSection_1,proj)
 ```
 
 for visualization, it is nice if we can remove the part of the cross-section that is above the topography.
-We can do that with the `BelowSurface` routine which returns a Boolean to indicate whether points are below or above the surface
+We can do that with the `belowSurface` routine which returns a Boolean to indicate whether points are below or above the surface
 
 ```julia
-below = BelowSurface(CrossSection_1_cart, TopoGeology_cart)
+below = belowSurface(CrossSection_1_cart, TopoGeology_cart)
 ```
 
 We can add that to the cross-section with:
@@ -240,14 +240,14 @@ Phases = zeros(Int8,size(ComputationalGrid.x)) #Define rock types
 Set everything below the topography to 1
 
 ```julia
-id = BelowSurface(ComputationalGrid, GeologyTopo_comp_surf)
+id = belowSurface(ComputationalGrid, GeologyTopo_comp_surf)
 Phases[id] .= 1
 ```
 
 The basement is set to 2
 
 ```julia
-id = BelowSurface(ComputationalGrid, Basement_comp_surf)
+id = belowSurface(ComputationalGrid, Basement_comp_surf)
 Phases[id] .= 2
 ```
 

docs/src/man/Tutorial_LaPalma.md

@@ -120,7 +120,7 @@ Phases = zeros(Int64,size(Grid_3D.x))
 Points that are below the surface are set to one:
 
 ```julia
-Below = BelowSurface(Grid_3D, Topo_model);
+Below = belowSurface(Grid_3D, Topo_model);
 Phases[Below] .= 1
 ```
 

docs/src/man/Tutorial_MohoTopo_Spada.md

@@ -0,0 +1,152 @@
+```@meta
+EditURL = "../../../tutorials/Tutorial_MohoTopo_Spada.jl"
+```
+
+# Moho Topography
+
+## Goal
+This explains how to load the Moho topography for Italy and the Alps and create a paraview file.
+The data comes from the following publication:
+    Spada, M., Bianchi, I., Kissling, E., Agostinetti, N.P., Wiemer, S., 2013. Combining controlled-source seismology and receiver function information to derive 3-D Moho topography for Italy. Geophysical Journal International 194, 1050–1068. doi:10.1093/gji/ggt148
+
+## 1. Download data
+The data is available as digital dataset on the Researchgate page of Prof. Edi Kissling
+[https://www.researchgate.net/publication/322682919\_Moho\_Map\_Data-WesternAlps-SpadaETAL2013](https://www.researchgate.net/publication/322682919_Moho_Map_Data-WesternAlps-SpadaETAL2013)
+
+We have also uploaded it here:
+[https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/](https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/)
+
+The full data set actually includes 3 different Moho's (Europe, Adria, Tyrrhenia-Corsica). To simplify matters, we have split the full file into 3 separate ascii files and uploaded it.
+
+We start with loading the necessary packages
+
+```julia
+using DelimitedFiles, GeophysicalModelGenerator
+```
+
+Please download the files
+- `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt`
+- `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt`
+- `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho3.txt`.
+This can be done using `download_data`:
+
+```julia
+download_data("https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/files/?p=%2FMoho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt&dl=1","Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt")
+download_data("https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/files/?p=%2FMoho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt&dl=1","Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt")
+download_data("https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/files/?p=%2FMoho_Map_Data-WesternAlps-SpadaETAL2013_Moho3.txt&dl=1","Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho3.txt")
+```
+
+## 2. Read data into Julia
+The data sets start at line 39. We read this into julia as:
+
+```julia
+data            = readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt",' ',Float64,'\n', skipstart=38,header=false)
+lon, lat, depth = data[:,1], data[:,2], -data[:,3];
+nothing #hide
+```
+
+Note that depth is made negative.
+
+## 3. Reformat the data
+Next, let's check if the data is spaced in a regular manner in Lon/Lat direction.
+For that, we plot it using the Plots package (you may have to install that first on your machine).
+
+```julia
+using Plots
+scatter(lon,lat,marker_z=depth, ylabel="latitude",xlabel="longitude",markersize=2.5, c = :roma)
+```
+
+![DataPoints](../assets/img/Tutorial_MohoSpada_LonLat.png)
+What we can see nicely here is that the data is reasonably regular but also that there are obviously locations where no data is define.
+
+The easiest way to transfer this to Paraview is to simply save this as 3D data points:
+
+```julia
+using GeophysicalModelGenerator
+data_Moho1 = GeoData(lon,lat,depth,(MohoDepth=depth*km,))
+```
+
+```julia
+GeoData
+  size  : (12355,)
+  lon   ϵ [ 4.00026 - 11.99991]
+  lat   ϵ [ 42.51778 - 48.99544]
+  depth ϵ [ -57.46 km - -21.34 km]
+  fields: (:MohoDepth,)
+```
+
+```julia
+Write_Paraview(data_Moho1, "Spada_Moho_Europe", PointsData=true)
+```
+
+And we can do the same with the other two Moho's:
+
+```julia
+data = readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt",' ',Float64,'\n', skipstart=38,header=false);
+lon, lat, depth        = data[:,1], data[:,2], -data[:,3];
+data_Moho2 = GeoData(lon,lat,depth,(MohoDepth=depth*km,))
+Write_Paraview(data_Moho2, "Spada_Moho_Adria", PointsData=true)
+data =readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho3.txt",' ',Float64,'\n', skipstart=38,header=false);
+lon, lat, depth        = data[:,1], data[:,2], -data[:,3];
+data_Moho3 = GeoData(lon,lat,depth,(MohoDepth=depth*km,))
+Write_Paraview(data_Moho3, "Spada_Moho_Tyrrhenia", PointsData=true)
+```
+
+If we plot this in paraview, it looks like this:
+![DataPoints_PV](../assets/img/Tutorial_MohoSpada_LonLat_Paraview.png)
+
+## 4. Fitting a mesh through the data
+So obviously, the Moho is discontinuous between these three Mohos. Often, it looks nicer if we fit a regular surface through these data points. To do this we first combine the data points of the 3 surfaces into one set of points
+
+```julia
+lon   = [data_Moho1.lon.val;   data_Moho2.lon.val;   data_Moho3.lon.val];
+lat   = [data_Moho1.lat.val;   data_Moho2.lat.val;   data_Moho3.lat.val];
+depth = [data_Moho1.depth.val; data_Moho2.depth.val; data_Moho3.depth.val];
+data_Moho_combined = GeoData(lon, lat, depth, (MohoDepth=depth*km,))
+```
+
+Next, we define a regular lon/lat grid
+
+```julia
+Lon,Lat,Depth = LonLatDepthGrid(4.1:0.1:11.9,42.5:.1:49,-30km);
+nothing #hide
+```
+
+And we will use a nearest neighbor interpolation method to fit a surface through the data. This has the advantage that it will take the discontinuities into account. We will use the package [NearestNeighbors.jl](https://github.com/KristofferC/NearestNeighbors.jl) for this, which you may have to install first
+
+```julia
+using NearestNeighbors
+kdtree = KDTree([lon'; lat'])
+idxs, dists = knn(kdtree, [Lon[:]'; Lat[:]'], 1, true)
+```
+
+`idxs` contains the indices of the closest points to the grid in (Lon,Lat). Next
+
+```julia
+Depth = zeros(size(Lon));
+for i=1:length(idxs)
+    Depth[i] = depth[idxs[i]][1]
+end
+```
+
+Now, we can create a `GeoData` structure with the regular surface and save it to paraview:
+
+```julia
+data_Moho = GeoData(Lon, Lat, Depth, (MohoDepth=Depth,))
+Write_Paraview(data_Moho, "Spada_Moho_combined")
+```
+
+The result is shown here, where the previous points are colored white and are a bit smaller. Obviously, the datasets coincide well.
+![DataPoints_Moho_surface](../assets/img/Tutorial_MohoSpada_Surface_Paraview.png)
+
+## 5. Julia script
+The full julia script that does it all is given [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/tree/main/tutorials/Tutorial_MohoTopo_Spada.jl). You need to be in the same directory as in the data file, after which you can run it in julia with
+
+```julia
+include("Tutorial_MohoTopo_Spada.jl")
+```
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+

docs/src/man/Tutorial_Votemaps.md

@@ -2,7 +2,7 @@
 EditURL = "../../../tutorials/Tutorial_Votemaps.jl"
 ```
 
-# 12 - Votemaps
+# Votemaps
 
 ## Aim
 In this tutorial, your will learn how to create Votemaps that compare different tomographic models and look for similarities between different models.