Merge pull request #68 from ChristianSchuler/add_Fault_Density_Tutorial

Add fault density tutorial; introduce CountMap function

docs/make.jl

@@ -96,7 +96,8 @@ makedocs(;
             "12 - VoteMaps" =>  "man/Tutorial_Votemaps.md",
             "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"
+            "15 - Create movies" =>  "man/tutorial_time_Seismicity.md",
+            "16 - Fault Density Map" =>  "man/tutorial_Fault_Map.md"
         ],
         "User Guide" => Any[
             "Installation" =>  "man/installation.md",

docs/src/man/tools.md

@@ -20,4 +20,5 @@ GeophysicalModelGenerator.Convert2CartData
 GeophysicalModelGenerator.ProjectCartData
 GeophysicalModelGenerator.DrapeOnTopo
 GeophysicalModelGenerator.LithostaticPressure!
+GeophysicalModelGenerator.CountMap
 ```

docs/src/man/tutorial_Fault_Map.md

@@ -0,0 +1,107 @@
+# 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).*
+

src/event_counts.jl

@@ -1,6 +1,6 @@
 using NearestNeighbors
 
-export PointData2NearestGrid
+export PointData2NearestGrid, CountMap
 
 
 """
@@ -114,3 +114,73 @@ function PointData2NearestGrid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)
     return count
 end
 
+"""
+    DatasetCountMap = CountMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64)
+
+Takes a 2D GeoData struct and counts entries of `field` per predefined control area. `field` should only consist of 1.0s and 0.0s. The control area is defined by `steplon` and `steplat`.
+`steplon` is the number of control areas in longitude direction and `steplat` the number of control areas in latitude direction.
+The counts per control area are normalized by the highest count.
+
+```julia
+julia> Data_Faults         = GeoData(Lon3D,Lat3D,Faults,(Faults=Faults,))
+GeoData 
+    size      : (375, 208, 1)
+    lon       ϵ [ -9.932408319802885 : 34.93985125012068]
+    lat       ϵ [ 35.086096468211394 : 59.919210145128545]
+    depth     ϵ [ 0.0 : 1.0]
+    fields    : (:Faults,)
+
+julia> steplon  = 125
+julia> steplat  = 70
+julia> countmap = CountMap(Data_Faults,"Faults",steplon,steplat)
+
+GeoData 
+    size      : (124, 69, 1)
+    lon       ϵ [ -9.751471789279 : 34.75891471959677]
+    lat       ϵ [ 35.26604656731949 : 59.73926004602028]
+    depth     ϵ [ 0.0 : 1.0]
+    fields    : (:CountMap,)
+```julia
+
+"""
+function CountMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64)
+
+    lon      = unique(DataSet.lon.val)
+    lat      = unique(DataSet.lat.val)
+
+    # create new lon/lat arrays which hold the boundaries of the control areas
+    lonstep  = LinRange(lon[1],lon[end],stepslon)
+    latstep  = LinRange(lat[1],lat[end],stepslat)
+    dlon     = abs(lonstep[2]-lonstep[1])
+    dlat     = abs(latstep[2]-latstep[1])
+    loncen   = lonstep[1]+dlon/2:dlon:lonstep[end]-dlon/2
+    latcen   = latstep[1]+dlat/2:dlat:latstep[end]-dlat/2
+    countmap = zeros(length(loncen),length(latcen))
+
+    expr =   Meta.parse(field)
+    if !haskey(DataSet.fields,expr[1])
+        error("The GeoData set does not have the field: $(expr[1])")
+    end
+
+    # count the ones in every control area
+    for i in eachindex(loncen)
+        for j in eachindex(latcen)
+            indi          = findall((lon .>= lonstep[i]) .& (lon .<= lonstep[i+1]))
+            indj          = findall((lat .>= latstep[j]) .& (lat .<= latstep[j+1]))
+            dataint       = DataSet.fields[expr[1]][indi,indj,1]
+            count         = sum(dataint)
+            countmap[i,j] = count
+        end
+    end
+
+    # normalize count in every control area
+    maxcount = maximum(countmap)
+    countmap = countmap ./ maxcount
+
+    # create new GeoData
+    Lon3D,Lat3D, Data = LonLatDepthGrid(loncen,latcen,0);
+    Data[:,:,1]       .= countmap
+    DatasetCountMap   = GeoData(Lon3D,Lat3D,Data,(CountMap=Data,))
+
+    return DatasetCountMap
+end

src/utils.jl

@@ -1,7 +1,7 @@
 # few utils that are useful
 
 export meshgrid, CrossSection, CrossSectionVolume, CrossSectionSurface, CrossSectionPoints, ExtractSubvolume, SubtractHorizontalMean
-export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap
+export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap, CountMap
 export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane
 export RotateTranslateScale
 export DrapeOnTopo, LithostaticPressure!

test/test_event_counts.jl

@@ -6,6 +6,11 @@ Grid_cart = CartData(XYZGrid(-20:20,-20:.1:20,-30:30))
 Lon,Lat,Depth = LonLatDepthGrid(-20:20,-20:.1:20,-30:30);
 Grid_geo = GeoData(Lon,Lat,Depth,(;Depth))
 
+# create 2D GeoData struct
+Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,0);
+CM              =   zeros(size(Depth)); CM[1:5,1:5] .= 1.0
+Data_set2D      =   GeoData(Lon,Lat,Depth,(Count=CM,))  
+
 using StableRNGs
 
 rng = StableRNG(123)
@@ -41,3 +46,9 @@ Grid_Count = PointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], Grid_geo; radius_fac
 # Test in case the EQ data is also specified as GeoData
 Grid_Count = PointData2NearestGrid(EQ_geo, Grid_geo; radius_factor=2)
 @test extrema(Grid_Count.fields.Count) == (0, 85)
+
+# Test CountMap
+Data_CountMap = CountMap(Data_set2D,"Count",5,5)
+@test Data_CountMap.fields.CountMap[1,1] == 1.0
+@test Data_CountMap.fields.CountMap[2,2] == 0.4444444444444444
+@test Data_CountMap.fields.CountMap[4,4] == 0.0

tutorials/Fault_Density.jl

@@ -0,0 +1,66 @@
+# # 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)