Update utils.jl

[mthielma]

Added a preselection of points to cross_section_points to reduce the number of e.g. earthquakes that have to be processed when a global dataset is used. (#160)

[github-actions]

Your PR requires formatting changes to meet the project's style guidelines.
Please consider running Runic (git runic main) to apply these changes.

Click here to view the suggested changes.

diff --git a/src/utils.jl b/src/utils.jl
index 2158df2..aba995f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -444,7 +444,7 @@ function cross_section_points(P::GeoData; Depth_level = nothing, Lat_level = not
     end
 
     if !isnothing(Lon_level)   # vertical slice @ given longitude
-        
+
         # to define the projection point, only choose events close to the desired profile
         p_Point = ProjectionPoint(Lat = sum(P.lat.val) / length(P.lat.val), Lon = Lon_level) # define the projection point (lat/lon) as the latitude and the mean of the longitudes of the data
         P_UTM = convert2UTMzone(P, p_Point) # convert to UTM
@@ -474,22 +474,22 @@ function cross_section_points(P::GeoData; Depth_level = nothing, Lat_level = not
         # convert P to UTM Data
         # to avoid projection issues, reduce the given point data set to a set that is located abound the desired profile
         # to choose the subset, we use a box around the profile with the profile width being added to the start and end points of the profile
-        # approximate formula: 
+        # approximate formula:
         # Latitude: 1 deg = 110.574 km --> 1 km = 1/110.574 deg
-        # Longitude: 1 deg = 111.320*cos(latitude) km --> 1km = 1/ 111.320 / cos(latitude) 
-        lat_add = 1.1 * ustrip(uconvert(u"km", section_width))/110.574 # multiply with 1.1 to make sure the box is large enough
-        lat_start = minimum([Start[2],End[2]]) - lat_add;
-        lat_end   = maximum([Start[2],End[2]])   + lat_add;
-        
-        lon_add = 1.1*ustrip(uconvert(u"km", section_width))/111.3209
+        # Longitude: 1 deg = 111.320*cos(latitude) km --> 1km = 1/ 111.320 / cos(latitude)
+        lat_add = 1.1 * ustrip(uconvert(u"km", section_width)) / 110.574 # multiply with 1.1 to make sure the box is large enough
+        lat_start = minimum([Start[2], End[2]]) - lat_add
+        lat_end = maximum([Start[2], End[2]]) + lat_add
+
+        lon_add = 1.1 * ustrip(uconvert(u"km", section_width)) / 111.3209
 
-        lon_start = minimum([Start[1],End[1]]) - lon_add/cos(deg2rad(minimum([Start[1],End[1]])))
-        lon_end   = maximum([Start[1],End[1]]) + lon_add/cos(deg2rad(maximum([Start[1],End[1]])))
+        lon_start = minimum([Start[1], End[1]]) - lon_add / cos(deg2rad(minimum([Start[1], End[1]])))
+        lon_end = maximum([Start[1], End[1]]) + lon_add / cos(deg2rad(maximum([Start[1], End[1]])))
 
-        ind = findall( (lon_start .< P.lon.val .< lon_end) .& (lat_start .< P.lat.val .< lat_end))
+        ind = findall((lon_start .< P.lon.val .< lon_end) .& (lat_start .< P.lat.val .< lat_end))
 
-        # now create a GeoData structure that only contains the subset, fields Magnitude and depth are hardcoded, more 
-        P_sub = GeoData(P.lon.val[ind], P.lat.val[ind], P.depth.val[ind], (Magnitude=P.fields.Magnitude[ind],Depth=P.fields.Depth[ind]))
+        # now create a GeoData structure that only contains the subset, fields Magnitude and depth are hardcoded, more
+        P_sub = GeoData(P.lon.val[ind], P.lat.val[ind], P.depth.val[ind], (Magnitude = P.fields.Magnitude[ind], Depth = P.fields.Depth[ind]))
         P_UTM = convert2UTMzone(P_sub, p_Point) # convert to UTM
 
         # create a GeoData set containing the points that create the profile plane (we need three points to uniquely define that plane)

[mthielma]

Tests are failing right now as the number of points after creating a profile are less than anticipated. I'm looking into that to see what is happening.

[albert-de-montserrat]

Suggested change

	ind = findall( (lon_start .< P.lon.val .< lon_end) .& (lat_start .< P.lat.val .< lat_end))
	ind = @. (lon_start < P.lon.val < lon_end) & (lat_start < P.lat.val < lat_end)

I think you can get away without findall here