diff --git a/AUTHORS.md b/AUTHORS.md index 4df1c401..cd9edfda 100644 --- a/AUTHORS.md +++ b/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 diff --git a/docs/make.jl b/docs/make.jl index ebec7750..d7452cf0 100644 --- a/docs/make.jl +++ b/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,7 +97,7 @@ 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[ @@ -105,6 +105,7 @@ makedocs(; "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", diff --git a/docs/src/man/Tutorial_FaultDensity.md b/docs/src/man/Tutorial_FaultDensity.md new file mode 100644 index 00000000..3a54b5a7 --- /dev/null +++ b/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).* + diff --git a/docs/src/man/Tutorial_Jura.md b/docs/src/man/Tutorial_Jura.md index e45ff9a9..25eee019 100644 --- a/docs/src/man/Tutorial_Jura.md +++ b/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 ``` diff --git a/docs/src/man/Tutorial_LaPalma.md b/docs/src/man/Tutorial_LaPalma.md index 0451f7a2..e7e37f31 100644 --- a/docs/src/man/Tutorial_LaPalma.md +++ b/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 ``` diff --git a/docs/src/man/Tutorial_MohoTopo_Spada.md b/docs/src/man/Tutorial_MohoTopo_Spada.md new file mode 100644 index 00000000..20fe1af2 --- /dev/null +++ b/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).* + diff --git a/docs/src/man/Tutorial_Votemaps.md b/docs/src/man/Tutorial_Votemaps.md index ce05f3b1..258bf675 100644 --- a/docs/src/man/Tutorial_Votemaps.md +++ b/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. diff --git a/docs/src/man/projection.md b/docs/src/man/projection.md index 4b777735..57ca5172 100644 --- a/docs/src/man/projection.md +++ b/docs/src/man/projection.md @@ -101,11 +101,10 @@ So this interpolates the topographic data from the `GeoData` to the orthogonal c You can do similar projections with full 3D data sets or pointwise data. -#### 3. List of functions +#### 3. List of relevant functions - -# ```@docs -# GeophysicalModelGenerator.Convert2CartData -# GeophysicalModelGenerator.ProjectCartData -# GeophysicalModelGenerator.Convert2UTMzone -# ``` +```@docs +GeophysicalModelGenerator.Convert2CartData +GeophysicalModelGenerator.ProjectCartData +GeophysicalModelGenerator.Convert2UTMzone +``` diff --git a/docs/src/man/surfaces.md b/docs/src/man/surfaces.md new file mode 100644 index 00000000..ac2f0bde --- /dev/null +++ b/docs/src/man/surfaces.md @@ -0,0 +1,15 @@ +# Surfaces + +We have a number of functions to deal with horizontal surfaces, which are defined as `GeoData` or `CartData` structures that have 3 dimensional `flat` coordinate array. Flat implies that the resolution is `n` by `m` by `1`, which is required to correctly display them in paraview. +The most straightforward surface is the topography (which you can obtain with `ImportTopo`). +Surfaces can be added and subtracted. + +```@docs +drape_on_topo +fit_surface_to_points +aboveSurface +belowSurface +interpolateDataOnSurface +is_surface +remove_NaN_Surface! +``` diff --git a/docs/src/man/tools.md b/docs/src/man/tools.md index b6791937..3744e620 100644 --- a/docs/src/man/tools.md +++ b/docs/src/man/tools.md @@ -8,9 +8,9 @@ GeophysicalModelGenerator.ExtractSubvolume GeophysicalModelGenerator.InterpolateDataFields GeophysicalModelGenerator.VoteMap GeophysicalModelGenerator.SubtractHorizontalMean -GeophysicalModelGenerator.AboveSurface -GeophysicalModelGenerator.BelowSurface -GeophysicalModelGenerator.InterpolateDataOnSurface +GeophysicalModelGenerator.aboveSurface +GeophysicalModelGenerator.belowSurface +GeophysicalModelGenerator.interpolateDataOnSurface GeophysicalModelGenerator.InterpolateTopographyOnPlane GeophysicalModelGenerator.ParseColumns_CSV_File GeophysicalModelGenerator.RotateTranslateScale! @@ -18,7 +18,7 @@ GeophysicalModelGenerator.PointData2NearestGrid GeophysicalModelGenerator.Convert2UTMzone GeophysicalModelGenerator.Convert2CartData GeophysicalModelGenerator.ProjectCartData -GeophysicalModelGenerator.DrapeOnTopo +GeophysicalModelGenerator.drape_on_topo GeophysicalModelGenerator.LithostaticPressure! GeophysicalModelGenerator.CountMap ``` diff --git a/docs/src/man/tutorial_Coastlines.md b/docs/src/man/tutorial_Coastlines.md index 675e9551..39c39ad2 100644 --- a/docs/src/man/tutorial_Coastlines.md +++ b/docs/src/man/tutorial_Coastlines.md @@ -1,4 +1,4 @@ -# 5 - Add coastlines +# Add coastlines ## Goal diff --git a/docs/src/man/tutorial_Fault_Map.md b/docs/src/man/tutorial_Fault_Map.md deleted file mode 100644 index f0d22c17..00000000 --- a/docs/src/man/tutorial_Fault_Map.md +++ /dev/null @@ -1,107 +0,0 @@ -# Fault Density Map - -## Aim - -````julia -#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 with the CountMap function -```` - -## Load Data - -Load packages - -````julia -using GeophysicalModelGenerator, Shapefile, Plots, Rasters, GeoDatasets, Interpolations -```` - -Data from "Active Faults of Eurasia Database AFEAD v2022" DOI:10.13140/RG.2.2.25509.58084 - -````julia -File = "AFEAD_v2022/AFEAD_v2022/AFEAD_v2022.shp" -```` - -Load data using Shapefile - -````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,)) -```` - -#### 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 -countmap = CountMap(Data_Faults,"Faults",steplon,steplat) -```` - -Plot the density map with coastlines - -````julia -lon = unique(countmap.lon.val) -lat = unique(countmap.lat.val) -coastlinesEurope = itp[lon,lat] -coastlinesEurope = map(y -> y > 1 ? 1 : y, coastlinesEurope) -heatmap(lon,lat,coastlinesEurope',colormap=cgrad(:gray1,rev=true),alpha=1.0); -heatmap!(lon,lat,countmap.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).* - diff --git a/docs/src/man/tutorial_GMT_Topography.md b/docs/src/man/tutorial_GMT_Topography.md index 7e0514ce..d233fb38 100644 --- a/docs/src/man/tutorial_GMT_Topography.md +++ b/docs/src/man/tutorial_GMT_Topography.md @@ -1,4 +1,4 @@ -# 4 - Extract topographic data from GMT.jl +# Extract topographic data from GMT.jl ## Goal diff --git a/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md b/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md index 48535876..c8abcaa1 100644 --- a/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md +++ b/docs/src/man/tutorial_GMT_Topography_GeologicalMap.md @@ -1,4 +1,4 @@ -# 8 - Extract ETOPO1 topographic data using GMT.jl and drape a geological map on top of the topography (given as raster graphics) +# Extract ETOPO1 topographic data using GMT.jl and drape a geological map on top of the topography (given as raster graphics) ## Goal diff --git a/docs/src/man/tutorial_ISC_data.md b/docs/src/man/tutorial_ISC_data.md index 07f6e9c6..8d43dad9 100644 --- a/docs/src/man/tutorial_ISC_data.md +++ b/docs/src/man/tutorial_ISC_data.md @@ -1,4 +1,4 @@ -# 10 - Plot ISC earthquake data +# Plot ISC earthquake data ## Goal This explains how to load earthquake data obtained from the ISC catalogue. diff --git a/docs/src/man/tutorial_MohoTopo.md b/docs/src/man/tutorial_MohoTopo.md deleted file mode 100644 index 16b607b2..00000000 --- a/docs/src/man/tutorial_MohoTopo.md +++ /dev/null @@ -1,115 +0,0 @@ -# 3 - Moho topography - -## Goal -This explains how to load the Moho topography for Italy and the Alps and create a paraview file - -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](https://doi.org/10.1093/gji/ggt148) - - -## Steps -#### 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. - -Please download the files `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt`, `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt` and `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 -julia> using DelimitedFiles -julia> data =readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt",' ',Float64,'\n', skipstart=38,header=false) -julia> lon, lat, depth = data[:,1], data[:,2], -data[:,3]; -``` -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 -julia> using Plots -julia> 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 -julia> using GeophysicalModelGenerator -julia> data_Moho1 = GeoData(lon,lat,depth,(MohoDepth=depth*km,)) -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 -julia> data =readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt",' ',Float64,'\n', skipstart=38,header=false); -julia> lon, lat, depth = data[:,1], data[:,2], -data[:,3]; -julia> data_Moho2 = GeoData(lon,lat,depth,(MohoDepth=depth*km,)) -julia> Write_Paraview(data_Moho2, "Spada_Moho_Adria", PointsData=true) -julia> data =readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho3.txt",' ',Float64,'\n', skipstart=38,header=false); -julia> lon, lat, depth = data[:,1], data[:,2], -data[:,3]; -julia> data_Moho3 = GeoData(lon,lat,depth,(MohoDepth=depth*km,)) -julia> 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) - - -#### 3.1 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]; -julia> lat = [data_Moho1.lat.val; data_Moho2.lat.val; data_Moho3.lat.val]; -julia> depth = [data_Moho1.depth.val; data_Moho2.depth.val; data_Moho3.depth.val]; -julia> data_Moho_combined = GeoData(lon, lat, depth, (MohoDepth=depth*km,)) -``` - -Next, we define a regular lon/lat grid -```julia -julia> Lon,Lat,Depth = LonLatDepthGrid(4.1:0.1:11.9,42.5:.1:49,-30km); -``` - -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 -julia> using NearestNeighbors -julia> kdtree = KDTree([lon'; lat']) -julia> idxs, dists = knn(kdtree, [Lon[:]'; Lat[:]'], 1, true) -``` -`idxs` contains the indices of the closest points to the grid in (Lon,Lat). Next -```julia -julia> Depth = zeros(size(Lon))*km; -julia> 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 -julia> data_Moho = GeoData(Lon, Lat, Depth, (MohoDepth=Depth,)) -julia> 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) - -#### 4. Julia script - -The full julia script that does it all is given [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/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 -julia> include("MohoTopo_Spada.jl") -``` diff --git a/docs/src/man/tutorial_Screenshot_To_Paraview.md b/docs/src/man/tutorial_Screenshot_To_Paraview.md index 5375dcd9..42dbd7b4 100644 --- a/docs/src/man/tutorial_Screenshot_To_Paraview.md +++ b/docs/src/man/tutorial_Screenshot_To_Paraview.md @@ -1,4 +1,4 @@ -# 6 - Import profiles/maps from published papers +# Import profiles/maps from published papers ## Goal Ideally, all data should be available in digital format, after which you could use the tools described in the other tutorial to transform them into `GeoData` and export them to VTK. @@ -7,7 +7,7 @@ Yet, the reality is different and often data is not (yet) available, or papers a For that reason, `GeophysicalModelGenerator` has tools that allow you to transfer a screenshot from any published paper into `GeoData/Paraview` and see it in 3D at the correct geographic location. This can be done for vertical profiles and for mapviews, which gives you a quick and easy way to see those papers in a new (3D) light. Here, we explain how. -- [6 - Import profiles/maps from published papers](#6---import-profilesmaps-from-published-papers) +- [Import profiles/maps from published papers](#import-profilesmaps-from-published-papers) - [Goal](#goal) - [General procedure](#general-procedure) - [1. Download data and crop images](#1-download-data-and-crop-images) diff --git a/docs/src/man/tutorial_UTM.md b/docs/src/man/tutorial_UTM.md index f6fd1753..e13ef20d 100644 --- a/docs/src/man/tutorial_UTM.md +++ b/docs/src/man/tutorial_UTM.md @@ -1,4 +1,4 @@ -# 11 - Read in UTM data +# Read in UTM data ## Goal diff --git a/docs/src/man/tutorial_load3DSeismicData.md b/docs/src/man/tutorial_load3DSeismicData.md index 92d9991b..48531e9e 100644 --- a/docs/src/man/tutorial_load3DSeismicData.md +++ b/docs/src/man/tutorial_load3DSeismicData.md @@ -1,4 +1,4 @@ -# 1 - 3D tomography model in CSV formation +# 3D tomography model in CSV formation ## Goal This explains how to load a 3D P-wave model and plot it in Paraview as a 3D volumetric data set. It also shows how you can create horizontal or vertical cross-sections through the data in a straightforward manner and how you can extract subsets of the data; @@ -109,7 +109,7 @@ julia> unique(lat[ind]) #### 3.1 Reshape data and save to paraview Next, we reshape the vectors with lon/lat/depth data into 3D matrixes: -``` +```julia julia> resolution = (length(unique(lon)), length(unique(lat)), length(unique(depth))) (121, 94, 101) julia> Lon = reshape(lon, resolution); @@ -124,13 +124,12 @@ julia> iz=findall(x -> x==-101.0, Depth[1,1,:]) 1-element Vector{Int64}: 91 julia> data=dVp_perc_3D[:,:,iz]; -julia> heatmap(unique(lon), unique(lat),data[:,:,1]', c=:roma,title="$(Depth[1,1,iz]) km") +heatmap(unique(lon), unique(lat),data[:,:,1]', c=:roma,title="$(Depth[1,1,iz]) km") ``` ![DataPoints](../assets/img/Tutorial_Zhao_LatLon_data_101km.png) - So this looks good. -Next we create a paraview file: +Next, we create a paraview file: ```julia julia> using GeophysicalModelGenerator julia> Data_set = GeoData(Lon,Lat,Depth,(dVp_Percentage=dVp_perc_3D,)) @@ -146,25 +145,21 @@ julia> Write_Paraview(Data_set, "Zhao_etal_2016_dVp_percentage") ``` #### 4. Plotting data in Paraview -In paraview you should open the `*.vts` file, and press `Apply` (left menu) after doing that. Once you did that you can select `dVp_Percentage` and `Surface` (see red ellipses below)/. +In paraview, you should open the `*.vts` file, and press `Apply` (left menu) after doing that. Once you did that you can select `dVp_Percentage` and `Surface` (see red ellipses below). In paraview you can open the file and visualize it. - ![Paraview_1](../assets/img/Tutorial_Zhao_Paraview_1.png) - This visualisation employs the default colormap, which is not particularly good. -You can change that by importing the roma colormap (using the link described earlier). For this, open the colormap editor and click the one with the heart on the right hand side. Next, import roma and select it. - +You can change that by importing the `roma` colormap (using the link described earlier). For this, open the colormap editor and click the one with the heart on the right hand side. Next, import `roma` and select it. ![Paraview_2](../assets/img/Tutorial_Zhao_Paraview_2.png) In order to change the colorrange select the button in the red ellipse and change the lower/upper bound. ![Paraview_3](../assets/img/Tutorial_Zhao_Paraview_3.png) -If you want to create a horizontal cross-section @ 200 km depth, you need to select the `Slice` tool, select `Sphere` as a clip type, set the center to `[0,0,0]` and set the radius to `6171` (=radius earth - 200 km). +If you want to create a horizontal cross-section at 200 km depth, you need to select the `Slice` tool, select `Sphere` as a clip type, set the center to `[0,0,0]` and set the radius to `6171` (=radius earth - 200 km). ![Paraview_4](../assets/img/Tutorial_Zhao_Paraview_4.png) - After pushing `Apply`, you'll see this: ![Paraview_5](../assets/img/Tutorial_Zhao_Paraview_5.png) @@ -173,9 +168,9 @@ If you want to plot iso-surfaces (e.g. at -3%), you can use the `Clip` option ag ![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. #### 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. There is a simple way to achieve this using the `CrossSection` function. @@ -210,7 +205,7 @@ As you see, this cross-section is not taken at exactly 10 degrees longitude. Tha If you wish to interpolate the data, specify `Interpolate=true`: ```julia -julia> Data_cross = CrossSection(Data_set, Lon_level=10, Interpolate=true) +julia> Data_cross = CrossSection(Data_set, Lon_level=10, Interpolate=true) GeoData size : (1, 100, 100) lon ϵ [ 10.0 : 10.0] @@ -234,7 +229,6 @@ julia> Write_Paraview(Data_cross, "Zhao_CrossSection_diagonal") ``` Here an image that shows the resulting cross-sections: - ![Paraview_7](../assets/img/Tutorial_Zhao_Paraview_7.png) #### 6. Extract a (3D) subset of the data @@ -265,26 +259,22 @@ GeoData fields: (:dVp_Percentage,) julia> Write_Paraview(Data_subset, "Zhao_Subset_interp") ``` -## Load and save data to disk -It would be useful to save the 3D data set we just created to disk, such that we can easily load it again at a later stage and create cross-sections etc, or compare it with other models. -It is quite easy to do so with the [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) package: +#### 7. Load and save data to disk +It would be useful to save the 3D data set we just created to disk, such that we can easily load it again at a later stage and create cross-sections etc, or compare it with other models. This can be done with the `save_GMG` function: ```julia -julia> using JLD2 -julia> jldsave("Zhao_Pwave.jld2"; Data_set) +julia> save_GMG("Zhao_Pwave", Data_set) ``` If you, at a later stage, want to load this file again do it as follows: ```julia -julia> using JLD2, GeophysicalModelGenerator -julia> Data_set_Zhao2016_Vp = load_object("Zhao_Pwave.jld2") +julia> Data_set_Zhao2016_Vp = load_GMG("Zhao_Pwave") ``` -## Julia script - +#### 8. Julia script The full julia script that does it all is given [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorial/Alps_VpModel_Zhao_etal_JGR2016.jl). You need to be in the same directory as in the data file, after which you can run it in julia with ```julia -julia> include("Alps_VpModel_Zhao_etal_JGR2016.jl") +include("Alps_VpModel_Zhao_etal_JGR2016.jl") ``` diff --git a/docs/src/man/tutorial_loadirregular3DSeismicData.md b/docs/src/man/tutorial_loadirregular3DSeismicData.md index fd0d21ec..3a5575e1 100644 --- a/docs/src/man/tutorial_loadirregular3DSeismicData.md +++ b/docs/src/man/tutorial_loadirregular3DSeismicData.md @@ -1,4 +1,4 @@ -# 7 - 3D tomography model that is re-interpolated on a regular grid +# 3D tomography model that is re-interpolated on a regular grid ## Goal This explains how to load a 3D seismic data set that is given in CSV format (comma separated ASCII), and plot it in paraview. The example is a shear-wave velocity model of the Alpine-Mediterranean region, described in: diff --git a/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md b/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md index 6ab8952e..cd5c67d4 100644 --- a/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md +++ b/docs/src/man/tutorial_loadregular3DSeismicData_netCDF.md @@ -1,4 +1,4 @@ -# 2 - 3D tomography model that is given as a netCDF file +# 3D tomography model that is given as a netCDF file ## Goal This explains how to load a 3D seismic data set that is given in netCDF format, and plot it in paraview. The example is a shear-wave velocity model of the Alpine-Mediterranean region, described in: @@ -12,12 +12,12 @@ The data is can be downloaded from [https://ds.iris.edu/files/products/emc/emc-f #### 2. Read data into Julia The main data-file, `El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc`, is given as netCDF file. To read in data of this type, it is necessary to load an appropriate package. Here, we will use the [https://github.com/JuliaGeo/NetCDF.jl](NetCDF.jl) package. Download and install the package with: - ```julia-repl + ```julia julia> using Pkg julia> Pkg.add("NetCDF") ``` First, let us have a look at the contents of this file (assuming that you are in the same directory where the file is located): - ```julia-repl + ```julia julia> using NetCDF julia> filename = ("El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc") julia> ncinfo("El-Sharkawy-etal-G3.2020-MeRE2020-Mediterranean-0.0.nc") @@ -92,7 +92,7 @@ latitude FLOAT latitude longitude FLOAT longitude Vs FLOAT longitude latitude depth ``` -Here we can see that there are four variables in this file, three of them (depth,latitude, longitude) having a single dimension and the fourth one (Vs) having dimensions of the three previous variables. The three one-dimensional vectors therefore denote a regular grid of coordinates defining the locations where Vs is stored. +Here, we can see that there are four variables in this file, three of them (`depth`,`latitude`, `longitude`) having a single dimension and the fourth one (`Vs`) having dimensions of the three previous variables. The three one-dimensional vectors therefore denote a regular grid of coordinates defining the locations where `Vs` is stored. To load this data, we can now simply use the command *ncread*: ```julia-repl julia> lat = ncread(filename,"latitude") @@ -104,7 +104,7 @@ julia> depth = -1 .* depth Note that we multiplied depth with -1. This is necessary to make depth to be negative, as that is what `GeophysicalModelGenerator.jl` expects. #### 3. Reformat the coordinate data -In the netCDF file, coordinates are given as 1D vectors denoting the location of nodes in a regular grid. However, `GeophysicalModelGenerator.jl` expects true 3D data, where each data point is assigned a latitude,longitude, depth and the respective property (here: Vs). To generate this full regular 3D grid, do the following: +In the netCDF file, coordinates are given as 1D vectors denoting the location of nodes in a regular grid. However, `GeophysicalModelGenerator.jl` expects true 3D data, where each data point is assigned a latitude,longitude, depth and the respective property (here: `Vs`). To generate this full regular 3D grid, do the following: ```julia-repl julia> using GeophysicalModelGenerator Lon3D,Lat3D,Depth3D = LonLatDepthGrid(lon, lat, depth); @@ -112,7 +112,7 @@ Lon3D,Lat3D,Depth3D = LonLatDepthGrid(lon, lat, depth); #### 4. Generate Paraview file Once the 3D coordinate matrix has been generated, producing a Paraview file is done with the following command ```julia -julia> Data_set = GeoData(Lon3D,Lat3D,Depth3D,(Vs_km_s=Vs_3D,)) +julia> Data_set = GeoData(Lon3D,Lat3D,Depth3D,(Vs_km_s=Vs_3D,)) GeoData size : (100, 100, 301) lon ϵ [ 29.0 - 51.0] @@ -137,7 +137,7 @@ If you want to clip the data set @ 200 km depth, you need to select the `Clip` t #### 6. Julia script -The full julia script that does it all is given [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorial/MeRe_ElSharkawy.jl). You need to be in the same directory as in the data file, after which you can run it in julia with +The full julia script that does it all is given [here](https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl/blob/main/tutorials/MeRe_ElSharkawy.jl). You need to be in the same directory as in the data file, after which you can run it in julia with ```julia julia> include("MeRe_ElSharkawy.jl") ``` diff --git a/docs/src/man/tutorial_local_Flegrei.md b/docs/src/man/tutorial_local_Flegrei.md index b1303bc3..cadda377 100644 --- a/docs/src/man/tutorial_local_Flegrei.md +++ b/docs/src/man/tutorial_local_Flegrei.md @@ -1,4 +1,4 @@ -# 13 - Campi Flegrei Volcano Tutorial using cartesian coordinates +# Campi Flegrei Volcano tutorial using cartesian coordinates ## Goal diff --git a/docs/src/man/tutorial_time_Seismicity.md b/docs/src/man/tutorial_time_Seismicity.md index dadc6a2e..14ffaa8e 100644 --- a/docs/src/man/tutorial_time_Seismicity.md +++ b/docs/src/man/tutorial_time_Seismicity.md @@ -1,4 +1,4 @@ -# 15 - How to create a movie that shows the temporal evolution of seismicity +# How to create a movie that shows the temporal evolution of seismicity ## Goal This tutorial creates a movie of the spatial variations in seismicity using the earthquakes previously visualized at Campi Flegrei caldera. We visualize it against the travel-time model and tomography: diff --git a/docs/src/man/visualise.md b/docs/src/man/visualise.md index 84755629..c1947ee8 100644 --- a/docs/src/man/visualise.md +++ b/docs/src/man/visualise.md @@ -64,4 +64,4 @@ julia> Visualise(Data_Cart, Topography=Topo_Cart) ```@docs Visualise -``` +``` \ No newline at end of file diff --git a/src/GeophysicalModelGenerator.jl b/src/GeophysicalModelGenerator.jl index 96e5cf6e..542c1dd3 100644 --- a/src/GeophysicalModelGenerator.jl +++ b/src/GeophysicalModelGenerator.jl @@ -33,6 +33,7 @@ export load include("data_types.jl") include("data_import.jl") include("coord_conversion.jl") +include("nearest_points.jl") include("utils.jl") include("Paraview_output.jl") include("transformation.jl") @@ -43,6 +44,7 @@ include("stl.jl") include("ProfileProcessing.jl") include("IO.jl") include("event_counts.jl") +include("surface_functions.jl") # Add optional routines (only activated when the packages are loaded) diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl index bad3ec91..92a93520 100644 --- a/src/LaMEM_io.jl +++ b/src/LaMEM_io.jl @@ -69,21 +69,21 @@ Creates a `CartData` struct from a LaMEM grid and from fields stored on that gri CartData(Grid::LaMEM_grid, fields::NamedTuple) = CartData(Grid.X, Grid.Y, Grid.Z, fields) """ - Below = BelowSurface(Data_LaMEM::LaMEM_grid, DataSurface_Cart::CartData) + Below = belowSurface(Data_LaMEM::LaMEM_grid, DataSurface_Cart::CartData) Determines if points within the 3D `LaMEM_grid` structure are below the Cartesian surface DataSurface_Cart """ -function BelowSurface(Grid::LaMEM_grid, DataSurface_Cart::CartData) - return AboveSurface(CartData(Grid,(Z=Grid.Z,)), DataSurface_Cart; above=false) +function belowSurface(Grid::LaMEM_grid, DataSurface_Cart::CartData) + return aboveSurface(CartData(Grid,(Z=Grid.Z,)), DataSurface_Cart; above=false) end """ - Above = AboveSurface(Data_LaMEM::LaMEM_grid, DataSurface_Cart::CartData) + Above = aboveSurface(Data_LaMEM::LaMEM_grid, DataSurface_Cart::CartData) Determines if points within the 3D `LaMEM_grid` structure are above the Cartesian surface DataSurface_Cart """ -function AboveSurface(Grid::LaMEM_grid, DataSurface_Cart::CartData) - return AboveSurface(CartData(Grid,(Z=Grid.Z,)), DataSurface_Cart; above=true) +function aboveSurface(Grid::LaMEM_grid, DataSurface_Cart::CartData) + return aboveSurface(CartData(Grid,(Z=Grid.Z,)), DataSurface_Cart; above=true) end diff --git a/src/data_types.jl b/src/data_types.jl index f783f291..768e9e13 100644 --- a/src/data_types.jl +++ b/src/data_types.jl @@ -1,7 +1,7 @@ # This is data_types.jl # contains type definitions to be used in GeophysicalModelGenerator -import Base: show +import Base: show, size export GeoData, ParaviewData, UTMData, CartData, LonLatDepthGrid, XYZGrid, Velocity_SphericalToCartesian!, @@ -216,6 +216,7 @@ struct GeoData <: AbstractGeneralGrid end end +size(d::GeoData) = size(d.lon.val) # Print an overview of the Geodata struct: function Base.show(io::IO, d::GeoData) @@ -247,6 +248,23 @@ function Base.show(io::IO, d::GeoData) end +""" + GeoData(lld::Tuple{Array,Array,Array}) + +This creates a `GeoData` struct if you have a Tuple with 3D coordinates as input. +# Example +```julia +julia> data = GeoData(LonLatDepthGrid(-10:10,-5:5,0)) +GeoData + size : (21, 11, 1) + lon ϵ [ -10.0 : 10.0] + lat ϵ [ -5.0 : 5.0] + depth ϵ [ 0.0 : 0.0] + fields : (:Z,) +``` +""" +GeoData(lld::Tuple) = GeoData(lld[1],lld[2],lld[3],(Z=lld[3],)) + """ ParaviewData(x::GeoUnit, y::GeoUnit, z::GeoUnit, values::NamedTuple) @@ -265,6 +283,7 @@ mutable struct ParaviewData <: AbstractGeneralGrid z :: GeoUnit fields :: NamedTuple end +size(d::ParaviewData) = size(d.x.val) # Print an overview of the ParaviewData struct: function Base.show(io::IO, d::ParaviewData) @@ -464,6 +483,7 @@ struct UTMData <: AbstractGeneralGrid end end +size(d::UTMData) = size(d.EW.val) # Print an overview of the UTMData struct: function Base.show(io::IO, d::UTMData) @@ -748,6 +768,7 @@ struct CartData <: AbstractGeneralGrid end end +size(d::CartData) = size(d.x.val) # Print an overview of the UTMData struct: function Base.show(io::IO, d::CartData) @@ -1032,6 +1053,7 @@ struct CartGrid{FT, D} <: AbstractGeneralGrid coord1D :: NTuple{D,StepRangeLen{FT}} # Tuple with 1D vectors in all directions coord1D_cen :: NTuple{D,StepRangeLen{FT}} # Tuple with 1D vectors of center points in all directions end +size(d::CartGrid) = d.N """ @@ -1149,8 +1171,6 @@ function CreateCartGrid(; end - - # view grid object function show(io::IO, g::CartGrid{FT, DIM}) where {FT, DIM} diff --git a/src/nearest_points.jl b/src/nearest_points.jl new file mode 100644 index 00000000..e0790d59 --- /dev/null +++ b/src/nearest_points.jl @@ -0,0 +1,68 @@ +# contains 1D,2D,3D nearest point routines +# mostly useful internally, but are exported as well (may be helpful in some cases) +using NearestNeighbors + +export nearest_point_indices + +# internal routine that does the computation + +""" + ind = nearest_point_indices(X::Array, X_pt::Vector) + +Returns the index of the nearest point in (`X_pt`) to (`X`) and returns the index +""" +function nearest_point_indices(X::Array, X_pt::Vector) + + # use nearest neighbour to interpolate data + coord = [X_pt';] + kdtree = KDTree(coord; leafsize = 10); + points = [vec(X)';]; + idx,_ = nn(kdtree, points); + + # transform to correct shape + ind = zeros(Int64,size(X)) + ind[:] = idx + + return ind +end + +""" + ind = nearest_point_indices(X::Array,Y::Array, X_pt::Vector, Y_pt::Vector) + +Returns the index of the nearest point in (`X_pt`,`Y_pt`) to (`X`,`Y`) and returns the index +""" +function nearest_point_indices(X::Array,Y::Array, X_pt::Vector, Y_pt::Vector) + + # use nearest neighbour to interpolate data + coord = [X_pt'; Y_pt']; + kdtree = KDTree(coord; leafsize = 10); + points = [vec(X)';vec(Y)']; + idx,_ = nn(kdtree, points); + + # transform to correct shape + ind = zeros(Int64,size(X)) + ind[:] = idx + + return ind +end + + +""" + ind = nearest_point_indices(X::Array,Y::Array,Z::Array, X_pt::Vector,Y_pt::Vector,Z_pt::Vector) + +Returns the index of the nearest point in (`X_pt`,`Y_pt`,`Z_pt`) to (`X`,`Y`,`Z`) and returns the index +""" +function nearest_point_indices(X::Array,Y::Array,Z::Array, X_pt::Vector,Y_pt::Vector,Z_pt::Vector) + + # use nearest neighbour to interpolate data + coord = [X_pt'; Y_pt'; Z_pt']; + kdtree = KDTree(coord; leafsize = 10); + points = [vec(X)';vec(Y)'; vec(Z)']; + idx,_ = nn(kdtree, points); + + # transform to correct shape + ind = zeros(Int64,size(X)) + ind[:] = idx + + return ind +end diff --git a/src/surface_functions.jl b/src/surface_functions.jl new file mode 100644 index 00000000..6d220c77 --- /dev/null +++ b/src/surface_functions.jl @@ -0,0 +1,397 @@ +# This contains a number of routines that deal with surfaces +export remove_NaN_Surface!, drape_on_topo, is_surface, fit_surface_to_points +export aboveSurface, belowSurface, interpolateDataOnSurface +import Base: +,- + +""" + issurf = is_surface(surf::AbstractGeneralGrid) + +Returns true if `surf` is a horizontal 3D surface. +""" +function is_surface(surf::AbstractGeneralGrid) + if size(surf)[3] == 1 + issurf = true + else + issurf = false + end + return issurf +end + +function +(a::_T, b::_T) where _T<:AbstractGeneralGrid + @assert size(a) == size(b) + return _addSurfaces(a,b) +end + +function -(a::_T, b::_T) where _T<:AbstractGeneralGrid + @assert size(a) == size(b) + return _subtractSurfaces(a,b) +end + +# Internal routines +_addSurfaces(a::_T, b::_T) where _T<:GeoData = GeoData(a.lon.val, a.lat.val, a.depth.val + b.depth.val, merge(a.fields,b.fields)) +_addSurfaces(a::_T, b::_T) where _T<:UTMData = UTMData(a.EW.val, a.NS.val, a.depth.val + b.depth.val, merge(a.fields,b.fields)) +_addSurfaces(a::_T, b::_T) where _T<:CartData = CartData(a.x.val, a.y.val, a.z.val + b.z.val, merge(a.fields,b.fields)) +_addSurfaces(a::_T, b::_T) where _T<:ParaviewData = ParaviewData(a.x.val, a.y.val, a.z.val + b.z.val, merge(a.fields,b.fields)) + +_subtractSurfaces(a::_T, b::_T) where _T<:GeoData = GeoData(a.lon.val, a.lat.val, a.depth.val - b.depth.val, merge(a.fields,b.fields)) +_subtractSurfaces(a::_T, b::_T) where _T<:UTMData = UTMData(a.EW.val, a.NS.val, a.depth.val - b.depth.val, merge(a.fields,b.fields)) +_subtractSurfaces(a::_T, b::_T) where _T<:CartData = CartData(a.x.val, a.y.val, a.z.val - b.z.val, merge(a.fields,b.fields)) +_subtractSurfaces(a::_T, b::_T) where _T<:ParaviewData = ParaviewData(a.x.val, a.y.val, a.z.val - b.z.val, merge(a.fields,b.fields)) + +""" + remove_NaN_Surface!(Z::Array,X::Array,Y::Array) + +Removes NaN's from a grid `Z` by taking the closest points as specified by `X` and `Y`. +""" +function remove_NaN_Surface!(Z,X,Y) + @assert size(Z) == size(X) == size(Y) + + # use nearest neighbour to interpolate data + id = findall(isnan.(Z) .== false) + id_NaN = findall(isnan.(Z)) + + coord = [X[id]'; Y[id]']; + kdtree = KDTree(coord; leafsize = 10); + + points = [X[id_NaN]'; Y[id_NaN]']; + idx,dist = nn(kdtree, points); + + Z[id_NaN] = Z[id[idx]] + + return nothing +end + + +""" + Topo = drape_on_topo(Topo::GeoData, Data::GeoData) + +This drapes fields of a data set `Data` on the topography `Topo` + +""" +function drape_on_topo(Topo::GeoData, Data::GeoData) + @assert is_surface(Topo) + @assert is_surface(Data) + + Lon,Lat,_ = LonLatDepthGrid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]); + + # use nearest neighbour to interpolate data + idx = nearest_point_indices(Lon,Lat, vec(Data.lon.val), vec(Data.lat.val) ); + + idx_out = findall( (Lon .< minimum(Data.lon.val)) .| (Lon .> maximum(Data.lon.val)) .| + (Lat .< minimum(Data.lat.val)) .| (Lat .> maximum(Data.lat.val)) ) + + fields_new = Topo.fields; + field_names = keys(Data.fields); + + for i = 1:length(Data.fields) + + if typeof(Data.fields[i]) <: Tuple + + # vector or anything that contains more than 1 field + data_tuple = Data.fields[i] # we have a tuple (likely a vector field), so we have to loop + + data_array = zeros(typeof(data_tuple[1][1]),size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3),length(Data.fields[i])); # create a 3D array that holds the 2D interpolated values + unit_array = zeros(size(data_array)); + + for j=1:length(data_tuple) + data_field = data_tuple[j]; + tmp = data_array[:,:,:,1]; + tmp = data_field[idx] + data_array[:,:,:,j] = tmp + end + + data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple + + # remove points outside domain + for j=1:length(data_tuple) + tmp = data_new[j]; + tmp[idx_out] .= NaN + data_array[:,:,:,j] = tmp + end + data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple + + else + # scalar field + data_new = zeros(typeof(Data.fields[i][1]), size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3)); + data_new = Data.fields[i][idx] # interpolate data field + + end + + # replace the one + new_field = NamedTuple{(field_names[i],)}((data_new,)) # Create a tuple with same name + fields_new = merge(fields_new, new_field); # replace the field in fields_new + + end + + Topo_new = GeoData(Topo.lon.val,Topo.lat.val,Topo.depth.val, fields_new) + + return Topo_new +end + + +""" + drape_on_topo(Topo::CartData, Data::CartData) + +Drapes Cartesian Data on topography +""" +function drape_on_topo(Topo::CartData, Data::CartData) + @assert is_surface(Topo) + @assert is_surface(Data) + + Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields ) + Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields ) + + Topo_new_lonlat = drape_on_topo(Topo_lonlat, Data_lonlat) + + Topo_new = CartData(Topo_new_lonlat.lon.val, Topo_new_lonlat.lat.val, Topo_new_lonlat.depth.val, Topo_new_lonlat.fields) + + return Topo_new +end + + +""" + surf_new = fit_surface_to_points(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector) + +This fits the `depth` values of the surface `surf` to the `depth` value of the closest-by-points in (`lon_pt`,`lat_pt`, `depth_pt`) + +""" +function fit_surface_to_points(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector) + @assert is_surface(surf) + + idx = nearest_point_indices(NumValue(surf.lon),NumValue(surf.lat), lon_pt, lat_pt); + depth = NumValue(surf.depth) + depth[idx] .= depth_pt[idx]; + + surf_new = surf + surf_new.depth .= depth + return surf_new +end + + +""" + surf_new = fit_surface_to_points(surf::CartData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector) + +This fits the `depth` values of the surface `surf` to the `depth` value of the closest-by-points in (`lon_pt`,`lat_pt`, `depth_pt`) + +""" +function fit_surface_to_points(surf::CartData, X_pt::Vector, Y_pt::Vector, Z_pt::Vector) + @assert is_surface(surf) + + idx = nearest_point_indices(NumValue(surf.x),NumValue(surf.y), X_pt[:], Y_pt[:]); + depth = NumValue(surf.z) + depth = Z_pt[idx] + + surf_new = deepcopy(surf) + surf_new.z.val .= depth + return surf_new +end + + + +""" + aboveSurface(Data::GeoData, DataSurface::GeoData; above=true) + +Returns a boolean array of size(Data.Lon), which is true for points that are above the surface DataSurface (or for points below if above=false). + +This can be used, for example, to mask points above/below the Moho in a volumetric dataset or in a profile. + +# Example +First we create a 3D data set and a 2D surface: +```julia +julia> Lon,Lat,Depth = LonLatDepthGrid(10:20,30:40,(-300:25:0)km); +julia> Data = Depth*2; +julia> Data_set3D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon)) +GeoData + size : (11, 11, 13) + lon ϵ [ 10.0 : 20.0] + lat ϵ [ 30.0 : 40.0] + depth ϵ [ -300.0 km : 0.0 km] + fields: (:Depthdata, :LonData) +julia> Lon,Lat,Depth = LonLatDepthGrid(10:20,30:40,-40km); +julia> Data_Moho = GeoData(Lon,Lat,Depth+Lon*km, (MohoDepth=Depth,)) + GeoData + size : (11, 11, 1) + lon ϵ [ 10.0 : 20.0] + lat ϵ [ 30.0 : 40.0] + depth ϵ [ -30.0 km : -20.0 km] + fields: (:MohoDepth,) +``` +Next, we intersect the surface with the data set: +```julia +julia> Above = aboveSurface(Data_set3D, Data_Moho); +``` +Now, `Above` is a boolean array that is true for points above the surface and false for points below and at the surface. + +""" +function aboveSurface(Data::GeoData, DataSurface::GeoData; above=true) + + if size(DataSurface.lon)[3]!=1 + error("It seems that DataSurface is not a surface") + end + + # Create interpolation object for surface + Lon_vec = DataSurface.lon.val[:,1,1]; + Lat_vec = DataSurface.lat.val[1,:,1]; + interpol = linear_interpolation((Lon_vec, Lat_vec), ustrip.(DataSurface.depth.val[:,:,1])); # create interpolation object + + DepthSurface = interpol.(Data.lon.val,Data.lat.val); + DepthSurface = DepthSurface*unit(DataSurface.depth.val[1]) + + if above + Above = Data.depth.val .> DepthSurface; + else + Above = Data.depth.val .< DepthSurface; + end + + return Above +end + +""" + Below = belowSurface(Data::GeoData, DataSurface::GeoData) + +Determines if points within the 3D `Data` structure are below the GeoData surface `DataSurface` +""" +function belowSurface(Data::GeoData, DataSurface::GeoData) + return aboveSurface(Data::GeoData, DataSurface::GeoData; above=false) +end + +""" + Above = aboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true) + +Determines if points within the 3D `Data_Cart` structure are above the Cartesian surface `DataSurface_Cart` +""" +function aboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true) + + Data = GeoData(ustrip.(Data_Cart.x.val), ustrip.(Data_Cart.y.val), ustrip.(Data_Cart.z.val), Data_Cart.fields) + DataSurface = GeoData(ustrip.(DataSurface_Cart.x.val),ustrip.(DataSurface_Cart.y.val), ustrip.(DataSurface_Cart.z.val), DataSurface_Cart.fields ) + + return Above = aboveSurface(Data, DataSurface; above=above) +end + +""" + Above = aboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true) + +Determines if points within the 3D `Data_Cart` structure are above the Cartesian surface `DataSurface_Cart` +""" +function aboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true) + + Data = GeoData(ustrip.(Data_Cart.x.val), ustrip.(Data_Cart.y.val), ustrip.(Data_Cart.z.val), Data_Cart.fields) + DataSurface = GeoData(ustrip.(DataSurface_Cart.x.val),ustrip.(DataSurface_Cart.y.val), ustrip.(DataSurface_Cart.z.val), DataSurface_Cart.fields ) + + return Above = aboveSurface(Data, DataSurface; above=above) +end + +""" + Above = aboveSurface(Grid::CartGrid, DataSurface_Cart::CartData; above=true) + +Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart` +""" +function aboveSurface(Grid::CartGrid, DataSurface_Cart::CartData; above=true) + + X,Y,Z = XYZGrid(Grid.coord1D...) + Data = CartData(Grid,(Z=Z,)) + + return aboveSurface(Data, DataSurface_Cart; above=above) +end + + +""" + Below = belowSurface(Grid::CartGrid, DataSurface_Cart::CartData) + + Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart` +""" +function belowSurface(Grid::CartGrid, DataSurface_Cart::CartData) + return aboveSurface(Grid, DataSurface_Cart; above=false) +end + + +""" + Below = belowSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData) + +Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart +""" +function belowSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData) + return aboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=false) +end + +""" + Below = belowSurface(Data_Cart::CartData, DataSurface_Cart::CartData) + +Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart +""" +function belowSurface(Data_Cart::CartData, DataSurface_Cart::CartData) + return aboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=false) +end + +""" + Surf_interp = interpolateDataOnSurface(V::ParaviewData, Surf::ParaviewData) + +Interpolates a 3D data set `V` on a surface defined by `Surf`. +# Example +```julia +julia> Data +ParaviewData + size : (33, 33, 33) + x ϵ [ -3.0 : 3.0] + y ϵ [ -2.0 : 2.0] + z ϵ [ -2.0 : 0.0] + fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res) +julia> surf +ParaviewData + size : (96, 96, 1) + x ϵ [ -2.9671875 : 3.2671875] + y ϵ [ -1.9791666666666667 : 1.9791666666666667] + z ϵ [ -1.5353766679763794 : -0.69925457239151] + fields: (:Depth,) +julia> Surf_interp = interpolateDataOnSurface(Data, surf) + ParaviewData + size : (96, 96, 1) + x ϵ [ -2.9671875 : 3.2671875] + y ϵ [ -1.9791666666666667 : 1.9791666666666667] + z ϵ [ -1.5353766679763794 : -0.69925457239151] + fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res) +``` +""" +function interpolateDataOnSurface(V::ParaviewData, Surf::ParaviewData) + + # Create GeoData structure: + V_geo = GeoData(V.x.val, V.y.val, V.z.val, V.fields) + V_geo.depth.val = ustrip(V_geo.depth.val); + + Surf_geo = GeoData(Surf.x.val, Surf.y.val, Surf.z.val, Surf.fields) + Surf_geo.depth.val = ustrip(Surf_geo.depth.val); + + Surf_interp_geo = interpolateDataOnSurface(V_geo, Surf_geo) + Surf_interp = ParaviewData(Surf_interp_geo.lon.val, Surf_interp_geo.lat.val, ustrip.(Surf_interp_geo.depth.val), Surf_interp_geo.fields) + + return Surf_interp + +end + +function interpolateDataOnSurface(V::CartData, Surf::CartData) + + # Create GeoData structure: + V_geo = GeoData(V.x.val, V.y.val, V.z.val, V.fields) + V_geo.depth.val = ustrip(V_geo.depth.val); + + Surf_geo = GeoData(Surf.x.val, Surf.y.val, Surf.z.val, Surf.fields) + Surf_geo.depth.val = ustrip(Surf_geo.depth.val); + + Surf_interp_geo = interpolateDataOnSurface(V_geo, Surf_geo) + Surf_interp = CartData(Surf_interp_geo.lon.val, Surf_interp_geo.lat.val, ustrip.(Surf_interp_geo.depth.val), Surf_interp_geo.fields) + + return Surf_interp + +end + +""" + Surf_interp = interpolateDataOnSurface(V::GeoData, Surf::GeoData) + +Interpolates a 3D data set `V` on a surface defined by `Surf` +""" +function interpolateDataOnSurface(V::GeoData, Surf::GeoData) + + Surf_interp = InterpolateDataFields(V, Surf.lon.val, Surf.lat.val, Surf.depth.val) + + return Surf_interp +end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 5561c1a5..4be325b9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,13 +1,12 @@ # few utils that are useful export meshgrid, CrossSection, CrossSectionVolume, CrossSectionSurface, CrossSectionPoints, ExtractSubvolume, SubtractHorizontalMean -export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap, CountMap -export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane +export ParseColumns_CSV_File, VoteMap, CountMap +export InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane export RotateTranslateScale -export DrapeOnTopo, LithostaticPressure! +export LithostaticPressure! export FlattenCrossSection export AddField, RemoveField -export SubtractSurfaces!, AddSurfaces!, remove_NaN_Surface! using NearestNeighbors @@ -1269,80 +1268,6 @@ function InterpolateDataFields2D_vecs(EW_vec, NS_vec, depth, fields_new, EW, NS) return depth_new, fields_new end - - -""" - Surf_interp = InterpolateDataOnSurface(V::ParaviewData, Surf::ParaviewData) - -Interpolates a 3D data set `V` on a surface defined by `Surf`. -# Example -```julia -julia> Data -ParaviewData - size : (33, 33, 33) - x ϵ [ -3.0 : 3.0] - y ϵ [ -2.0 : 2.0] - z ϵ [ -2.0 : 0.0] - fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res) -julia> surf -ParaviewData - size : (96, 96, 1) - x ϵ [ -2.9671875 : 3.2671875] - y ϵ [ -1.9791666666666667 : 1.9791666666666667] - z ϵ [ -1.5353766679763794 : -0.69925457239151] - fields: (:Depth,) -julia> Surf_interp = InterpolateDataOnSurface(Data, surf) - ParaviewData - size : (96, 96, 1) - x ϵ [ -2.9671875 : 3.2671875] - y ϵ [ -1.9791666666666667 : 1.9791666666666667] - z ϵ [ -1.5353766679763794 : -0.69925457239151] - fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res) -``` -""" -function InterpolateDataOnSurface(V::ParaviewData, Surf::ParaviewData) - - # Create GeoData structure: - V_geo = GeoData(V.x.val, V.y.val, V.z.val, V.fields) - V_geo.depth.val = ustrip(V_geo.depth.val); - - Surf_geo = GeoData(Surf.x.val, Surf.y.val, Surf.z.val, Surf.fields) - Surf_geo.depth.val = ustrip(Surf_geo.depth.val); - - Surf_interp_geo = InterpolateDataOnSurface(V_geo, Surf_geo) - Surf_interp = ParaviewData(Surf_interp_geo.lon.val, Surf_interp_geo.lat.val, ustrip.(Surf_interp_geo.depth.val), Surf_interp_geo.fields) - - return Surf_interp - -end - -function InterpolateDataOnSurface(V::CartData, Surf::CartData) - - # Create GeoData structure: - V_geo = GeoData(V.x.val, V.y.val, V.z.val, V.fields) - V_geo.depth.val = ustrip(V_geo.depth.val); - - Surf_geo = GeoData(Surf.x.val, Surf.y.val, Surf.z.val, Surf.fields) - Surf_geo.depth.val = ustrip(Surf_geo.depth.val); - - Surf_interp_geo = InterpolateDataOnSurface(V_geo, Surf_geo) - Surf_interp = CartData(Surf_interp_geo.lon.val, Surf_interp_geo.lat.val, ustrip.(Surf_interp_geo.depth.val), Surf_interp_geo.fields) - - return Surf_interp - -end - -""" - Surf_interp = InterpolateDataOnSurface(V::GeoData, Surf::GeoData) - -Interpolates a 3D data set `V` on a surface defined by `Surf` -""" -function InterpolateDataOnSurface(V::GeoData, Surf::GeoData) - - Surf_interp = InterpolateDataFields(V, Surf.lon.val, Surf.lat.val, Surf.depth.val) - - return Surf_interp -end # Extracts a sub-data set using indices function ExtractDataSets(V::AbstractGeneralGrid, iLon, iLat, iDepth) @@ -1467,10 +1392,6 @@ function SubtractHorizontalMean(V::AbstractArray{T, 2}; Percentage=false) where return V_sub end - - - - """ ParseColumns_CSV_File(data_file, num_columns) @@ -1509,143 +1430,6 @@ function ParseColumns_CSV_File(data_file, num_columns) return data end - -""" - AboveSurface(Data::GeoData, DataSurface::GeoData; above=true) - -Returns a boolean array of size(Data.Lon), which is true for points that are above the surface DataSurface (or for points below if above=false). - -This can be used, for example, to mask points above/below the Moho in a volumetric dataset or in a profile. - -# Example -First we create a 3D data set and a 2D surface: -```julia -julia> Lon,Lat,Depth = LonLatDepthGrid(10:20,30:40,(-300:25:0)km); -julia> Data = Depth*2; -julia> Data_set3D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon)) -GeoData - size : (11, 11, 13) - lon ϵ [ 10.0 : 20.0] - lat ϵ [ 30.0 : 40.0] - depth ϵ [ -300.0 km : 0.0 km] - fields: (:Depthdata, :LonData) -julia> Lon,Lat,Depth = LonLatDepthGrid(10:20,30:40,-40km); -julia> Data_Moho = GeoData(Lon,Lat,Depth+Lon*km, (MohoDepth=Depth,)) - GeoData - size : (11, 11, 1) - lon ϵ [ 10.0 : 20.0] - lat ϵ [ 30.0 : 40.0] - depth ϵ [ -30.0 km : -20.0 km] - fields: (:MohoDepth,) -``` -Next, we intersect the surface with the data set: -```julia -julia> Above = AboveSurface(Data_set3D, Data_Moho); -``` -Now, `Above` is a boolean array that is true for points above the surface and false for points below and at the surface. - -""" -function AboveSurface(Data::GeoData, DataSurface::GeoData; above=true) - - if size(DataSurface.lon)[3]!=1 - error("It seems that DataSurface is not a surface") - end - - # Create interpolation object for surface - Lon_vec = DataSurface.lon.val[:,1,1]; - Lat_vec = DataSurface.lat.val[1,:,1]; - interpol = linear_interpolation((Lon_vec, Lat_vec), ustrip.(DataSurface.depth.val[:,:,1])); # create interpolation object - - DepthSurface = interpol.(Data.lon.val,Data.lat.val); - DepthSurface = DepthSurface*unit(DataSurface.depth.val[1]) - - if above - Above = Data.depth.val .> DepthSurface; - else - Above = Data.depth.val .< DepthSurface; - end - - return Above -end - -""" - Below = BelowSurface(Data::GeoData, DataSurface::GeoData) - -Determines if points within the 3D `Data` structure are below the GeoData surface `DataSurface` -""" -function BelowSurface(Data::GeoData, DataSurface::GeoData) - return AboveSurface(Data::GeoData, DataSurface::GeoData; above=false) -end - -""" - Above = AboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true) - -Determines if points within the 3D `Data_Cart` structure are above the Cartesian surface `DataSurface_Cart` -""" -function AboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true) - - Data = GeoData(ustrip.(Data_Cart.x.val), ustrip.(Data_Cart.y.val), ustrip.(Data_Cart.z.val), Data_Cart.fields) - DataSurface = GeoData(ustrip.(DataSurface_Cart.x.val),ustrip.(DataSurface_Cart.y.val), ustrip.(DataSurface_Cart.z.val), DataSurface_Cart.fields ) - - return Above = AboveSurface(Data, DataSurface; above=above) -end - -""" - Above = AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true) - -Determines if points within the 3D `Data_Cart` structure are above the Cartesian surface `DataSurface_Cart` -""" -function AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true) - - Data = GeoData(ustrip.(Data_Cart.x.val), ustrip.(Data_Cart.y.val), ustrip.(Data_Cart.z.val), Data_Cart.fields) - DataSurface = GeoData(ustrip.(DataSurface_Cart.x.val),ustrip.(DataSurface_Cart.y.val), ustrip.(DataSurface_Cart.z.val), DataSurface_Cart.fields ) - - return Above = AboveSurface(Data, DataSurface; above=above) -end - -""" - Above = AboveSurface(Grid::CartGrid, DataSurface_Cart::CartData; above=true) - -Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart` -""" -function AboveSurface(Grid::CartGrid, DataSurface_Cart::CartData; above=true) - - X,Y,Z = XYZGrid(Grid.coord1D...) - Data = CartData(Grid,(Z=Z,)) - - return AboveSurface(Data, DataSurface_Cart; above=above) -end - - -""" - Below = BelowSurface(Grid::CartGrid, DataSurface_Cart::CartData) - - Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart` -""" -function BelowSurface(Grid::CartGrid, DataSurface_Cart::CartData) - return AboveSurface(Grid, DataSurface_Cart; above=false) -end - - -""" - Below = BelowSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData) - -Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart -""" -function BelowSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData) - return AboveSurface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=false) -end - -""" - Below = BelowSurface(Data_Cart::CartData, DataSurface_Cart::CartData) - -Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart -""" -function BelowSurface(Data_Cart::CartData, DataSurface_Cart::CartData) - return AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=false) -end - - """ VoteMap(DataSets::Vector{GeoData}, criteria::Vector{String}, dims=(50,50,50)) @@ -1841,96 +1625,6 @@ function RotateTranslateScale(Data::Union{ParaviewData, CartData}; Rotate::Numbe end -""" - Topo = DrapeOnTopo(Topo::GeoData, Data::GeoData) - -This drapes fields of a data set `Data` on the topography `Topo` - - -""" -function DrapeOnTopo(Topo::GeoData, Data::GeoData) - - - Lon,Lat,Depth = LonLatDepthGrid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]); - - # use nearest neighbour to interpolate data - coord = [vec(Data.lon.val)'; vec(Data.lat.val)']; - kdtree = KDTree(coord; leafsize = 10); - points = [vec(Lon)';vec(Lat)']; - idx,dist = nn(kdtree, points); - - - idx_out = findall( (Lon .< minimum(Data.lon.val)) .| (Lon .> maximum(Data.lon.val)) .| - (Lat .< minimum(Data.lat.val)) .| (Lat .> maximum(Data.lat.val)) ) - - fields_new = Topo.fields; - field_names = keys(Data.fields); - - for i = 1:length(Data.fields) - - if typeof(Data.fields[i]) <: Tuple - - # vector or anything that contains more than 1 field - data_tuple = Data.fields[i] # we have a tuple (likely a vector field), so we have to loop - - data_array = zeros(typeof(data_tuple[1][1]),size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3),length(Data.fields[i])); # create a 3D array that holds the 2D interpolated values - unit_array = zeros(size(data_array)); - - for j=1:length(data_tuple) - data_field = data_tuple[j]; - tmp = data_array[:,:,:,1]; - tmp = data_field[idx] - data_array[:,:,:,j] = tmp - end - - data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple - - # remove points outside domain - for j=1:length(data_tuple) - tmp = data_new[j]; - tmp[idx_out] .= NaN - data_array[:,:,:,j] = tmp - end - data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple - - else - - # scalar field - data_new = zeros(typeof(Data.fields[i][1]), size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3)); - data_new = Data.fields[i][idx] # interpolate data field - - end - - # replace the one - new_field = NamedTuple{(field_names[i],)}((data_new,)) # Create a tuple with same name - fields_new = merge(fields_new, new_field); # replace the field in fields_new - - end - - - Topo_new = GeoData(Topo.lon.val,Topo.lat.val,Topo.depth.val, fields_new) - - return Topo_new - -end - - -""" - DrapeOnTopo(Topo::CartData, Data::CartData) - -Drapes Cartesian Data on topography -""" -function DrapeOnTopo(Topo::CartData, Data::CartData) - Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields ) - Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields ) - - Topo_new_lonlat = DrapeOnTopo(Topo_lonlat, Data_lonlat) - - Topo_new = CartData(Topo_new_lonlat.lon.val, Topo_new_lonlat.lat.val, Topo_new_lonlat.depth.val, Topo_new_lonlat.fields) - - return Topo_new -end - """ LithostaticPressure!(Plithos::Array, Density::Array, dz::Number; g=9.81) @@ -1946,83 +1640,4 @@ function LithostaticPressure!(Plithos::Array{T,N}, Density::Array{T,N}, dz::Numb Plithos[:] = reverse!(cumsum(reverse!(Plithos),dims=N)) return nothing -end - - -""" - AddSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - -Adds `Surface2` to `Surface1`. The addition happens on the `Surface1.depth`; the fields remain unchanged - -""" -function AddSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - - Surface1.depth.val .= Surface1.depth.val + Surface2.depth.val - - return nothing -end - -""" - AddSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - -Adds `Surface2` to `Surface1`. The addition happens on the `Surface1.z`; the fields remain unchanged - -""" -function AddSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - - Surface1.z.val .= Surface1.z.val + Surface2.z.val - - return nothing -end - - -""" - SubtractSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - -Subtracts `Surface2` to `Surface1`. The addition happens on the `Surface1.depth`; the fields remain unchanged - -""" -function SubtractSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - - Surface1.depth.val .= Surface1.depth.val - Surface2.depth.val - - return nothing -end - - -""" - SubtractSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - -Subtracts `Surface2` to `Surface1`. The addition happens on the `Surface1.z`; the fields remain unchanged - -""" -function SubtractSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - - Surface1.z.val .= Surface1.z.val - Surface2.z.val - - return nothing -end - - - - -""" - remove_NaN_Surface!(Z,X,Y) - -Removes NaN's from a grid `Z` by taking the closest points as specified by `X` and `Y`. -""" -function remove_NaN_Surface!(Z,X,Y) - # use nearest neighbour to interpolate data - id = findall(isnan.(Z) .== false) - id_NaN = findall(isnan.(Z)) - - coord = [X[id]'; Y[id]']; - kdtree = KDTree(coord; leafsize = 10); - - points = [X[id_NaN]'; Y[id_NaN]']; - idx,dist = nn(kdtree, points); - - Z[id_NaN] = Z[id[idx]] - - return nothing -end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 90ec22bc..ee6d9c42 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,12 +13,18 @@ end @testset "Gravity model" begin include("test_voxel_gravity.jl") end +@testset "Nearest points" begin + include("test_nearest_points.jl") +end @testset "Utils" begin include("test_utils.jl") end @testset "Transformations" begin include("test_transformation.jl") end +@testset "Surfaces" begin + include("test_surfaces.jl") +end @testset "LaMEM" begin include("test_lamem.jl") diff --git a/test/test_nearest_points.jl b/test/test_nearest_points.jl new file mode 100644 index 00000000..18d294c7 --- /dev/null +++ b/test/test_nearest_points.jl @@ -0,0 +1,21 @@ +using Test, GeophysicalModelGenerator + +# 3D arrays +c_1D = CartData(XYZGrid(1:4,0,0)) +c_2D = CartData(XYZGrid(1:4,1:5,2)) +c_3D = CartData(XYZGrid(1:4,1:5,2:5)) + +# points +X_pt, Y_pt, Z_pt = XYZGrid(1:.05:5,0:.07:8,1:.4:5) + +# 1D test +id_1D = nearest_point_indices(NumValue(c_1D.x),X_pt[:]) +@test sum(id_1D) == 166336 + +# 2D test +id_2D = nearest_point_indices(NumValue(c_2D.x),NumValue(c_2D.y), X_pt[:], Y_pt[:]) +@test sum(id_2D) == 992141 + +# 3D test +id_3D = nearest_point_indices(NumValue(c_2D.x),NumValue(c_2D.y),NumValue(c_2D.z), X_pt[:],Y_pt[:],Z_pt[:]) +@test sum(id_3D) == 442556 \ No newline at end of file diff --git a/test/test_setup_geometry.jl b/test/test_setup_geometry.jl index a6477d1f..b7834692 100644 --- a/test/test_setup_geometry.jl +++ b/test/test_setup_geometry.jl @@ -76,22 +76,6 @@ Phases = Compute_Phase(Phases, Temp, Grid, LP); @test Phases[1,1,25] == 1 @test Phases[1,1,73] == 0 - -# test AboveSurface with the Grid object -Grid = CreateCartGrid(size=(10,20,30),x=(0.,10), y=(0.,10), z=(-10.,2.)) -@test Grid.Δ[2] ≈ 0.5263157894736842 - -Temp = ones(Float64, Grid.N...)*1350; -Phases = zeros(Int32, Grid.N...); - - -Topo_cart = CartData(XYZGrid(-1:.2:20,-12:.2:13,0)); -ind = AboveSurface(Grid, Topo_cart); -@test sum(ind[1,1,:]) == 5 - -ind = BelowSurface(Grid, Topo_cart); -@test sum(ind[1,1,:]) == 25 - # Create Grid & nondimensionalize it CharDim = GEO_units(); Grid = CreateCartGrid(size=(10,20,30),x=(0.0km,10km), y=(0.0km, 10km), z=(-10.0km, 2.0km), CharDim=CharDim) diff --git a/test/test_surfaces.jl b/test/test_surfaces.jl new file mode 100644 index 00000000..82446f32 --- /dev/null +++ b/test/test_surfaces.jl @@ -0,0 +1,84 @@ +using Test +# test various surface routines + +# Create surfaces +cartdata1 = CartData(XYZGrid(1:4,1:5,0)) +cartdata2 = CartData(XYZGrid(1:4,1:5,2)) +cartdata3 = CartData(XYZGrid(1:4,1:5,2:5)) +cartdata2 = AddField(cartdata2,"Z2",cartdata2.x.val) + +@test is_surface(cartdata1) +@test is_surface(cartdata2) +@test is_surface(cartdata3) == false + +geodata1 = GeoData(LonLatDepthGrid(1:4,1:5,0)) +geodata2 = GeoData(LonLatDepthGrid(1:4,1:5,2)) +geodata3 = GeoData(LonLatDepthGrid(1:4,1:5,2:5)) + +@test is_surface(geodata1) +@test is_surface(geodata2) +@test is_surface(geodata3) == false + +# Test add & subtraction of surfaces +cartdata4 = cartdata1 + cartdata2 +@test length(cartdata4.fields)==2 +@test cartdata4.z.val[2]==2.0 + +cartdata5 = cartdata1 - cartdata2 +@test length(cartdata5.fields)==2 +@test cartdata5.z.val[2]==-2.0 + +geodata4 = geodata1 + geodata2 +@test length(geodata4.fields)==1 +@test geodata4.depth.val[2]==2.0 + +geodata5 = geodata1 - geodata2 +@test length(geodata5.fields)==1 +@test geodata5.depth.val[2]==-2.0 + +# Test removing NaN; +Z = NumValue(cartdata5.z) +Z[2,2] = NaN; +remove_NaN_Surface!(Z,NumValue(cartdata5.x), NumValue(cartdata5.y)) +@test any(isnan.(Z))==false + +# Test draping values on topography +X,Y,Z = XYZGrid(1:.14:4,1:.02:5,0); +v = X.^2 .+ Y.^2; +values1 = CartData(X,Y,Z, (; v)) +values2 = CartData(X,Y,Z, (; colors=(v,v,v) )) + +cart_drape1 = drape_on_topo(cartdata2, values1) +@test sum(cart_drape1.fields.v) ≈ 366.02799999999996 + +cart_drape2 = drape_on_topo(cartdata2, values2) +@test cart_drape2.fields.colors[1][10] ≈ 12.9204 + +values1 = GeoData(X,Y,Z, (; v)) +values2 = GeoData(X,Y,Z, (; colors=(v,v,v) )) + +geo_drape1 = drape_on_topo(geodata2, values1) +@test sum(geo_drape1.fields.v) ≈ 366.02799999999996 + +geo_drape2 = drape_on_topo(geodata2, values2) +@test geo_drape2.fields.colors[1][10] ≈ 12.9204 + +# test fit_surface_to_points +cartdata2b = fit_surface_to_points(cartdata2, X[:], Y[:], v[:]) +@test sum(NumValue(cartdata2b.z)) ≈ 366.02799999999996 + + +#------------- +# test aboveSurface with the Grid object +Grid = CreateCartGrid(size=(10,20,30),x=(0.,10), y=(0.,10), z=(-10.,2.)) +@test Grid.Δ[2] ≈ 0.5263157894736842 + +Temp = ones(Float64, Grid.N...)*1350; +Phases = zeros(Int32, Grid.N...); + +Topo_cart = CartData(XYZGrid(-1:.2:20,-12:.2:13,0)); +ind = aboveSurface(Grid, Topo_cart); +@test sum(ind[1,1,:]) == 5 + +ind = belowSurface(Grid, Topo_cart); +@test sum(ind[1,1,:]) == 25 diff --git a/test/test_utils.jl b/test/test_utils.jl index e2f28e83..b0dd62cb 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -168,15 +168,15 @@ Data_Moho = GeoData(Lon,Lat,Depth,(MohoDepth=Depth,LonData=Lon,TestData= # Test intersecting a surface with 2D or 3D data sets -Above = AboveSurface(Data_set3D, Data_Moho); # 3D regular ordering +Above = aboveSurface(Data_set3D, Data_Moho); # 3D regular ordering @test Above[1,1,12]==true @test Above[1,1,11]==false -Above = AboveSurface(Data_set3D_reverse, Data_Moho); # 3D reverse depth ordering +Above = aboveSurface(Data_set3D_reverse, Data_Moho); # 3D reverse depth ordering @test Above[1,1,2]==true @test Above[1,1,3]==false -Above = AboveSurface(Data_sub_cross, Data_Moho); # 2D cross-section +Above = aboveSurface(Data_sub_cross, Data_Moho); # 2D cross-section @test Above[end]==true @test Above[1]==false diff --git a/tutorials/Fault_Density.jl b/tutorials/Fault_Density.jl deleted file mode 100644 index b25104e6..00000000 --- a/tutorials/Fault_Density.jl +++ /dev/null @@ -1,66 +0,0 @@ -# # Fault Density Map - -# ## Aim -#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 with the CountMap function - -# ## Load Data - -# Load packages - -using GeophysicalModelGenerator, Shapefile, Plots, Rasters, GeoDatasets, Interpolations - -# Data from "Active Faults of Eurasia Database AFEAD v2022" DOI:10.13140/RG.2.2.25509.58084 -File = "AFEAD_v2022/AFEAD_v2022/AFEAD_v2022.shp" - -# Load data using Shapefile - -table = Shapefile.Table(File) -geoms = Shapefile.shapes(table) -CONF = table.CONF - -# Raster the shapefile data -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 -lonC,latC,dataC = GeoDatasets.landseamask(;resolution='l',grid=10) - -# Interpolate to fault grid -itp = linear_interpolation((lonC, latC), dataC) -coastlines = itp[lon.val,lat.val] -coastlines = map(y -> y > 1 ? 1 : y, coastlines) - -# Plot the fault data -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 -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 -Lon3D,Lat3D, Faults = LonLatDepthGrid(Lon,Lat,0); -Faults[:,:,1] = data -Data_Faults = GeoData(Lon3D,Lat3D,Faults,(Faults=Faults,)) - -# #### 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. -steplon = 188 -steplat = 104 -countmap = CountMap(Data_Faults,"Faults",steplon,steplat) - -# Plot the density map with coastlines -lon = unique(countmap.lon.val) -lat = unique(countmap.lat.val) -coastlinesEurope = itp[lon,lat] -coastlinesEurope = map(y -> y > 1 ? 1 : y, coastlinesEurope) -heatmap(lon,lat,coastlinesEurope',colormap=cgrad(:gray1,rev=true),alpha=1.0); -heatmap!(lon,lat,countmap.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) diff --git a/tutorials/MohoTopo_Spada.jl b/tutorials/MohoTopo_Spada.jl deleted file mode 100644 index 09e51c5d..00000000 --- a/tutorials/MohoTopo_Spada.jl +++ /dev/null @@ -1,48 +0,0 @@ -# This shows how to read in moho topography, plot the data as points and fit a surface through the data - -# 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](https://doi.org/10.1093/gji/ggt148) -# -# We have uploaded it here: https://seafile.rlp.net/d/a50881f45aa34cdeb3c0/ -# -# Please download the files `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt`, `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt` and `Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho3.txt` - -using DelimitedFiles, GeophysicalModelGenerator - -# read European Moho -data = readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt",' ',Float64,'\n', skipstart=38,header=false) -lon, lat, depth = data[:,1], data[:,2], -data[:,3]; -data_Moho1 = GeoData(lon,lat,depth,(MohoDepth=depth*km,)) - -# Write as data points: -Write_Paraview(data_Moho1, "Spada_Moho_Europe", PointsData=true) - -# Do the same with the other Moho's: -data =readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho2.txt",' ',Float64,'\n', skipstart=38,header=false); -lon, lat, depth = data[:,1], data[:,2], -data[:,3]; -Tutorial_MohoSpada_LonLat_Paraview = 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) - -# Combine all data points into one data set: -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,)) - -# Define regular grid -Lon,Lat,Depth = LonLatDepthGrid(4.1:0.1:11.9,42.5:.1:49,-30km); - -# Interpolate using NN: -using NearestNeighbors -kdtree = KDTree([lon'; lat']) -idxs, dists = knn(kdtree, [Lon[:]'; Lat[:]'], 1, true) -Depth = zeros(size(Lon))*km; -for i=1:length(idxs) - Depth[i] = depth[idxs[i]][1] -end -data_Moho = GeoData(Lon, Lat, Depth, (MohoDepth=Depth,)) -Write_Paraview(data_Moho, "Spada_Moho_combined") - diff --git a/tutorials/Tutorial_FaultDensity.jl b/tutorials/Tutorial_FaultDensity.jl new file mode 100644 index 00000000..375ea669 --- /dev/null +++ b/tutorials/Tutorial_FaultDensity.jl @@ -0,0 +1,70 @@ +# # 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: +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: +File = "AFEAD_v2022/AFEAD_v2022.shp" + +# Load data using the `Shapefile` package: +table = Shapefile.Table(File) +geoms = Shapefile.shapes(table) +CONF = table.CONF + +# Raster the shapefile data +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`: +lonC,latC,dataC = GeoDatasets.landseamask(;resolution='l',grid=10) + +# Interpolate to fault grid +itp = linear_interpolation((lonC, latC), dataC) +coastlines = itp[lon.val,lat.val] +coastlines = map(y -> y > 1 ? 1 : y, coastlines) + +# Plot the fault data +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 +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 +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. +steplon = 188 +steplat = 104 +cntmap = CountMap(Data_Faults,"Faults",steplon,steplat) + +# Plot the density map with coastlines +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`: +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) + + +#src Literate.markdown("Tutorial_FaultDensity.jl","docs/src/man",keepcomments=true, execute=false, codefence = "```julia" => "```") diff --git a/tutorials/Tutorial_Jura.jl b/tutorials/Tutorial_Jura.jl index 0171991d..1120f321 100644 --- a/tutorials/Tutorial_Jura.jl +++ b/tutorials/Tutorial_Jura.jl @@ -1,4 +1,4 @@ -# # 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. @@ -27,7 +27,7 @@ upperright = [8.948117154811715, 47.781282316442606, 0.0] Geology = Screenshot_To_GeoData("SchoriM_Encl_01_Jura-map_A1.png", lowerleft, upperright, fieldname=:geology_colors) # name should have "colors" in it # You can "drape" this image on the topographic map with -TopoGeology = DrapeOnTopo(Topo, Geology) +TopoGeology = drape_on_topo(Topo, Geology) # ```julia # GeoData # size : (3721, 2641, 1) @@ -48,7 +48,7 @@ download_data("https://zenodo.org/records/10726801/files/BMes_Spline_longlat.tif # Now, import the GeoTIFF as: Basement = ImportGeoTIFF("BMes_Spline_longlat.tif", fieldname=:Basement, removeNaN_z=true) # 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. @@ -125,8 +125,8 @@ Basement_cart = ProjectCartData(TopoGeology_cart, Basement, proj) 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 -below = BelowSurface(CrossSection_1_cart, TopoGeology_cart) +# We can do that with the `belowSurface` routine which returns a Boolean to indicate whether points are below or above the surface +below = belowSurface(CrossSection_1_cart, TopoGeology_cart) # We can add that to the cross-section with: CrossSection_1_cart = AddField(CrossSection_1_cart,"rocks",Int64.(below)) @@ -164,11 +164,11 @@ Basement_comp_surf = InterpolateDataFields2D(Basement_cart_rot, Computatio Phases = zeros(Int8,size(ComputationalGrid.x)) #Define rock types # Set everything below the topography to 1 -id = BelowSurface(ComputationalGrid, GeologyTopo_comp_surf) +id = belowSurface(ComputationalGrid, GeologyTopo_comp_surf) Phases[id] .= 1 # The basement is set to 2 -id = BelowSurface(ComputationalGrid, Basement_comp_surf) +id = belowSurface(ComputationalGrid, Basement_comp_surf) Phases[id] .= 2 # Add to the computational grid: diff --git a/tutorials/Tutorial_LaPalma.jl b/tutorials/Tutorial_LaPalma.jl index 7c1d0bb7..c4f53ac1 100644 --- a/tutorials/Tutorial_LaPalma.jl +++ b/tutorials/Tutorial_LaPalma.jl @@ -71,7 +71,7 @@ Grid_3D =PointData2NearestGrid(EQ_cart, Grid_3D, radius_factor=3) Phases = zeros(Int64,size(Grid_3D.x)) # Points that are below the surface are set to one: -Below = BelowSurface(Grid_3D, Topo_model); +Below = belowSurface(Grid_3D, Topo_model); Phases[Below] .= 1 # Lets assume that the crust is 15 km thick diff --git a/tutorials/Tutorial_MohoTopo_Spada.jl b/tutorials/Tutorial_MohoTopo_Spada.jl new file mode 100644 index 00000000..acaea6e9 --- /dev/null +++ b/tutorials/Tutorial_MohoTopo_Spada.jl @@ -0,0 +1,100 @@ +# # 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 +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`: +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: +data = readdlm("Moho_Map_Data-WesternAlps-SpadaETAL2013_Moho1.txt",' ',Float64,'\n', skipstart=38,header=false) +lon, lat, depth = data[:,1], data[:,2], -data[:,3]; +# 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). +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: +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,) +# ``` +Write_Paraview(data_Moho1, "Spada_Moho_Europe", PointsData=true) + +# And we can do the same with the other two Moho's: +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 +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 +Lon,Lat,Depth = LonLatDepthGrid(4.1:0.1:11.9,42.5:.1:49,-30km); + +# 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 +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 +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: + +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 +include("Tutorial_MohoTopo_Spada.jl") diff --git a/tutorials/Tutorial_Votemaps.jl b/tutorials/Tutorial_Votemaps.jl index 1aebf5d6..597f26cc 100644 --- a/tutorials/Tutorial_Votemaps.jl +++ b/tutorials/Tutorial_Votemaps.jl @@ -1,4 +1,4 @@ -# # 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.