event_counts.jl

Author: boriskaus

docs/src/man/Tutorial_FaultDensity.md

@@ -82,12 +82,12 @@ 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.
+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)
+cntmap   = countMap(Data_Faults,"Faults",steplon,steplat)
 ```
 
 Plot the density map with coastlines
@@ -103,7 +103,7 @@ 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")
+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)

docs/src/man/Tutorial_LaPalma.md

@@ -108,7 +108,7 @@ Grid_3D = CartData(XYZGrid(-35:.3:30,-15:.25:45,-50:.5:5))
 Next we check how many earthquakes are around the grid points:
 
 ```julia
-Grid_3D =PointData2NearestGrid(EQ_cart, Grid_3D, radius_factor=3)
+Grid_3D =pointData2NearestGrid(EQ_cart, Grid_3D, radius_factor=3)
 ```
 
 And we can define an array with rock types:

docs/src/man/tools.md

@@ -14,11 +14,11 @@ interpolateDataOnSurface
 InterpolateTopographyOnPlane
 ParseColumns_CSV_File
 RotateTranslateScale!
-PointData2NearestGrid
+pointData2NearestGrid
 Convert2UTMzone
 Convert2CartData
 projectCartData
 drape_on_topo
 LithostaticPressure!
-CountMap
+countMap
 ```

src/event_counts.jl

@@ -1,10 +1,10 @@
 using NearestNeighbors
 
-export PointData2NearestGrid, CountMap
+export pointData2NearestGrid, countMap
 
 
 """
-    Grid_counts = PointData2NearestGrid(Point::CartData, Grid::CartData; radius_factor=1)
+    Grid_counts = pointData2NearestGrid(Point::CartData, Grid::CartData; radius_factor=1)
 
 Uses nearest neighbour interpolation to count how many points (given by `Point`) are in the vicinity of a 3D `Grid`. 
 The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
@@ -13,38 +13,38 @@ The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
 
 `Grid_counts` is `Grid` but with an additional field `Count` that has the number of hits
 """
-function PointData2NearestGrid(Point::CartData, Grid::CartData; radius_factor=1)
+function pointData2NearestGrid(Point::CartData, Grid::CartData; radius_factor=1)
 
     @assert length(size(Point.x)) == 1
 
     # call routine
-    Count = PointData2NearestGrid(NumValue(Point.x),NumValue(Point.y), NumValue(Point.z), NumValue(Grid.x),NumValue(Grid.y),NumValue(Grid.z); radius_factor=radius_factor)
+    Count = pointData2NearestGrid(NumValue(Point.x),NumValue(Point.y), NumValue(Point.z), NumValue(Grid.x),NumValue(Grid.y),NumValue(Grid.z); radius_factor=radius_factor)
 
     # return CartGrid with added field
     return  AddField(Grid,"Count",Count);
 end
 
 
 """
-    Grid_counts = PointData2NearestGrid(pt_x,pt_y,pt_z, Grid::CartData; radius_factor=1)
+    Grid_counts = pointData2NearestGrid(pt_x,pt_y,pt_z, Grid::CartData; radius_factor=1)
 
 Uses nearest neighbour interpolation to count how many points (given by `pt_x`,`pt_y`,`pt_z` coordinate vectors) are in the 
 vicinity of 3D `CartGrid` specified by `Grid`. The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
 
 `Grid_counts` is `Grid` but with an additional field `Count` that has the number of hits
 """
-function PointData2NearestGrid(pt_x,pt_y,pt_z, Grid::CartData; radius_factor=1)
+function pointData2NearestGrid(pt_x,pt_y,pt_z, Grid::CartData; radius_factor=1)
 
     # call routine
-    Count = PointData2NearestGrid(pt_x,pt_y,pt_z, NumValue(Grid.x),NumValue(Grid.y),NumValue(Grid.z); radius_factor=radius_factor)
+    Count = pointData2NearestGrid(pt_x,pt_y,pt_z, NumValue(Grid.x),NumValue(Grid.y),NumValue(Grid.z); radius_factor=radius_factor)
 
     # return CartGrid with added field
     return  AddField(Grid,"Count",Count);
 end
 
 
 """
-    Grid_counts = PointData2NearestGrid(Point::GeoData, Grid::GeoData; radius_factor=1)
+    Grid_counts = pointData2NearestGrid(Point::GeoData, Grid::GeoData; radius_factor=1)
 
 Uses nearest neighbour interpolation to count how many points (given by `Point`) are in the vicinity of a 3D `Grid`. 
 The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
@@ -53,44 +53,44 @@ The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
 
 `Grid_counts` is `Grid` but with an additional field `Count` that has the number of hits
 """
-function PointData2NearestGrid(Point::GeoData, Grid::GeoData; radius_factor=1)
+function pointData2NearestGrid(Point::GeoData, Grid::GeoData; radius_factor=1)
 
     @assert length(size(Point.lon)) == 1
 
     # call routine
-    Count = PointData2NearestGrid(NumValue(Point.lon),NumValue(Point.lat), NumValue(Point.depth), NumValue(Grid.lon),NumValue(Grid.lat),NumValue(Grid.depth); radius_factor=radius_factor)
+    Count = pointData2NearestGrid(NumValue(Point.lon),NumValue(Point.lat), NumValue(Point.depth), NumValue(Grid.lon),NumValue(Grid.lat),NumValue(Grid.depth); radius_factor=radius_factor)
 
     # return CartGrid with added field
     return  AddField(Grid,"Count",Count);
 end
 
 
 """
-    Grid_counts = PointData2NearestGrid(pt_x,pt_y,pt_z, Grid::GeoData; radius_factor=1)
+    Grid_counts = pointData2NearestGrid(pt_x,pt_y,pt_z, Grid::GeoData; radius_factor=1)
 
 Uses nearest neighbour interpolation to count how many points (given by `pt_x`,`pt_y`,`pt_z` coordinate vectors) are in the 
 vicinity of 3D `GeoData` specified by `Grid`. The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
 
 `Grid_counts` is `Grid` but with an additional field `Count` that has the number of hits
 """
-function PointData2NearestGrid(pt_x,pt_y,pt_z, Grid::GeoData; radius_factor=1)
+function pointData2NearestGrid(pt_x,pt_y,pt_z, Grid::GeoData; radius_factor=1)
 
     # call routine
-    Count = PointData2NearestGrid(pt_x,pt_y,pt_z, NumValue(Grid.lon),NumValue(Grid.lat),NumValue(Grid.depth); radius_factor=radius_factor)
+    Count = pointData2NearestGrid(pt_x,pt_y,pt_z, NumValue(Grid.lon),NumValue(Grid.lat),NumValue(Grid.depth); radius_factor=radius_factor)
 
     # return CartGrid with added field
     return  AddField(Grid,"Count",Count);
 end
 
 """
-    count = PointData2NearestGrid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)
+    count = pointData2NearestGrid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)
 
 This uses nearest neighbour interpolation to count how many points (given by `pt_x`,`pt_y`,`pt_z` coordinate vectors) are in the 
 vicinity of 3D grid point specified by `X`,`Y`,`Z` 3D coordinate arrays, with regular spacing `(Δx,Δy,Δz)`.
 The search radius is `R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)`
 
 """
-function PointData2NearestGrid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)
+function pointData2NearestGrid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)
 
     data = zeros(3,length(pt_x));
     data[1,:],data[2,:],data[3,:] = pt_x[:], pt_y[:], pt_z[:]
@@ -115,7 +115,7 @@ function PointData2NearestGrid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)
 end
 
 """
-    DatasetCountMap = CountMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64)
+    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.
@@ -132,18 +132,18 @@ GeoData
 
 julia> steplon  = 125
 julia> steplat  = 70
-julia> countmap = CountMap(Data_Faults,"Faults",steplon,steplat)
+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,)
+    fields    : (:countMap,)
 ```julia
 
 """
-function CountMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64)
+function countMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64)
 
     lon      = unique(DataSet.lon.val)
     lat      = unique(DataSet.lat.val)
@@ -155,7 +155,7 @@ function CountMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64
     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))
+    countMap = zeros(length(loncen),length(latcen))
 
     expr =   Meta.parse(field)
     if !haskey(DataSet.fields,expr[1])
@@ -169,18 +169,18 @@ function CountMap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64
             indj          = findall((lat .>= latstep[j]) .& (lat .<= latstep[j+1]))
             dataint       = DataSet.fields[expr[1]][indi,indj,1]
             count         = sum(dataint)
-            countmap[i,j] = count
+            countMap[i,j] = count
         end
     end
 
     # normalize count in every control area
-    maxcount = maximum(countmap)
-    countmap = countmap ./ maxcount
+    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,))
+    Data[:,:,1]       .= countMap
+    DatasetcountMap   = GeoData(Lon3D,Lat3D,Data,(countMap=Data,))
 
-    return DatasetCountMap
+    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, VoteMap, CountMap
+export ParseColumns_CSV_File, VoteMap, countMap
 export InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane
 export RotateTranslateScale
 export LithostaticPressure!

test/test_event_counts.jl

@@ -28,27 +28,27 @@ R = (sum(pt.^2,dims=2)).^(1/3)
 ind = findall(R.< 5);
 
 # test the basic routine
-counts = PointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], NumValue(Grid_cart.x),NumValue(Grid_cart.y), NumValue(Grid_cart.z); radius_factor=2)
+counts = pointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], NumValue(Grid_cart.x),NumValue(Grid_cart.y), NumValue(Grid_cart.z); radius_factor=2)
 @test extrema(counts) == (0, 85)
 
 # Test if the grid is on a CartData grid
-Grid_Count = PointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], Grid_cart; radius_factor=2)
+Grid_Count = pointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], Grid_cart; radius_factor=2)
 @test extrema(Grid_Count.fields.Count) == (0, 85)
 
 # Test in case the EQ data is also specified as CartData
-Grid_Count = PointData2NearestGrid(EQ_cart, Grid_cart; radius_factor=2)
+Grid_Count = pointData2NearestGrid(EQ_cart, Grid_cart; radius_factor=2)
 @test extrema(Grid_Count.fields.Count) == (0, 85)
 
 # Test if the grid is on a GeoData grid
-Grid_Count = PointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], Grid_geo; radius_factor=2)
+Grid_Count = pointData2NearestGrid(pt[:,1],pt[:,2],pt[:,3], Grid_geo; radius_factor=2)
 @test extrema(Grid_Count.fields.Count) == (0, 85)
 
 # Test in case the EQ data is also specified as GeoData
-Grid_Count = PointData2NearestGrid(EQ_geo, Grid_geo; radius_factor=2)
+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
+# 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/Tutorial_FaultDensity.jl

@@ -50,10 +50,10 @@ 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.
+# 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)
+cntmap   = countMap(Data_Faults,"Faults",steplon,steplat)
 
 # Plot the density map with coastlines
 lon = unique(cntmap.lon.val)
@@ -63,7 +63,7 @@ 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")
+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)
 
 

tutorials/Tutorial_LaPalma.jl

@@ -65,7 +65,7 @@ Write_Paraview(Topo_model,"Topo_model")
 Grid_3D = CartData(XYZGrid(-35:.3:30,-15:.25:45,-50:.5:5))
 
 # Next we check how many earthquakes are around the grid points:
-Grid_3D =PointData2NearestGrid(EQ_cart, Grid_3D, radius_factor=3)
+Grid_3D =pointData2NearestGrid(EQ_cart, Grid_3D, radius_factor=3)
 
 # And we can define an array with rock types:
 Phases = zeros(Int64,size(Grid_3D.x))