Merge branch 'main' into mt_tutorial

Project.toml

@@ -15,6 +15,7 @@ Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MeshIO = "7269a6da-0436-5bbc-96c2-40638cbb6118"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
@@ -59,4 +60,4 @@ GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test","GMT"]
+test = ["Test", "GMT"]

src/data_import.jl

@@ -4,7 +4,9 @@
 #
 # Author: Marcel Thielmann, 05/2021
 
-export Screenshot_To_GeoData, Screenshot_To_CartData, Screenshot_To_UTMData
+using LightXML
+
+export Screenshot_To_GeoData, Screenshot_To_CartData, Screenshot_To_UTMData, GetLonLatDepthMag_QuakeML
 
 # import CSV data using standard library functions
 # here we assume that the data is indeed comma separated and that comments are preceded with a "#"
@@ -309,7 +311,43 @@ function Screenshot_To_UTMData(filename::String, Corner_LowerLeft, Corner_UpperR
 
       # first create a GeoData struct
       Data_UTM = Screenshot_To_GeoData(filename, Corner_LowerLeft, Corner_UpperRight; Corner_LowerRight=Corner_LowerRight, Corner_UpperLeft=Corner_UpperLeft, Cartesian=false, UTM=true, UTMzone=UTMzone, isnorth=isnorth, fieldname=fieldname)
-
       return Data_UTM
+end
 
+"""
+    Data = GetLonLatDepthMag_QuakeML(filename::String)
+
+Extracts longitude, latitude, depth and magnitude from a QuakeML file that has been e.g. downloaded from ISC. The data is then returned in GeoData format.
+"""
+function GetLonLatDepthMag_QuakeML(filename::String)
+    # The QuakeML format consists of a tree with quite a lot of branches, so we have to traverse it to quite some extent to get the desired values
+    # using LightXML: extension???
+    xdoc = parse_file(filename); # parse the whole file
+    xroot =root(xdoc);
+    catalogues = get_elements_by_tagname(xroot,"eventParameters");
+    catalogue  = catalogues[1];
+    events     = get_elements_by_tagname(catalogue,"event"); # now those are all events
+    num_events = size(events,1);
+
+    # allocate, lat,lon,depth,magnitude
+    lon    = zeros(num_events,1);
+    lat    = zeros(num_events,1);
+    depth  = zeros(num_events,1);
+    mag    = zeros(num_events,1);
+
+    # now loop over the events and assign the respective values
+    for ievent = 1:num_events
+        tmp_event    = events[ievent];
+        origin = get_elements_by_tagname(events[ievent], "origin");
+        magnitude = get_elements_by_tagname(events[ievent], "magnitude");
+
+        # this is a bit dirty, if you find a better/cleaner way, be my guest...
+        lon[ievent]   = parse(Float64,string(collect(child_nodes(collect(child_elements(get_elements_by_tagname(origin[1], "longitude")[1]))[1]))[1]))
+        lat[ievent]   = parse(Float64,string(collect(child_nodes(collect(child_elements(get_elements_by_tagname(origin[1], "latitude")[1]))[1]))[1]))
+        depth[ievent] = parse(Float64,string(collect(child_nodes(collect(child_elements(get_elements_by_tagname(origin[1], "depth")[1]))[1]))[1]))
+        mag[ievent]   = parse(Float64,string(collect(child_nodes(get_elements_by_tagname(get_elements_by_tagname(magnitude[1],"mag")[1],"value")[1]))[1]));
+    end
+
+    Data_ISC = GeoData(lon,lat,-1*depth/1e3,(Magnitude=mag,Depth=-1*depth/1e3*km));
+    return Data_ISC
 end

src/utils.jl

@@ -1173,6 +1173,23 @@ function InterpolateDataFields2D(Original::CartData, New::CartData; Rotate=0.0,
     return CartData(New.x.val,New.y.val,Znew, fields_new)
 end
 
+"""    
+    Surf_interp = InterpolateDataFields2D(V::GeoData, x::AbstractRange, y::AbstractRange;  Lat::Number, Lon::Number)
+
+Interpolates a 3D data set `V` with a projection point `proj=(Lat, Lon)` on a plane defined by `x` and `y`, where `x` and `y` are uniformly spaced.
+Returns the 2D array `Surf_interp`. 
+"""
+function InterpolateDataFields2D(V::GeoData, x::AbstractRange, y::AbstractRange;  Lat=49.9929, Lon=8.2473)
+    # Default: Lat=49.9929, Lon=8.2473 => Mainz (center of universe)
+    proj = ProjectionPoint(; Lat = Lat, Lon = Lon)
+    return InterpolateDataFields2D(V::GeoData, proj, x, y)
+end
+
+function InterpolateDataFields2D(LonLat::GeoData, proj::ProjectionPoint, x::AbstractRange, y::AbstractRange)
+    cart_grid = CartData(XYZGrid(x, y, 0))
+    tproj = ProjectCartData(cart_grid, LonLat, proj)
+    return tproj.z.val[:, :, 1]
+end
 
 """
     InterpolateDataFields2D_vecs(x_vec, y_vec, depth, fields_new, X, Y)
@@ -1326,19 +1343,7 @@ function InterpolateDataOnSurface(V::GeoData, Surf::GeoData)
 
     return Surf_interp
 end
-
-"""    
-    Surf_interp = InterpolateTopographyOnPlane(V::GeoData, proj::ProjectionPoint, x::AbstractRange, y::AbstractRange)
-
-Interpolates a 3D data set `V` with a projection point `proj` on a plane defined by `x` and `y`, where are uniformly spaced.
-Returns the 2D array `Surf_interp`. 
-"""
-function InterpolateTopographyOnPlane(V::GeoData, proj::ProjectionPoint, x::AbstractRange, y::AbstractRange)
-    cart_grid = CartData(XYZGrid(x, y, 0))
-    tproj = ProjectCartData(cart_grid, V, proj)
-    return tproj.z.val[:, :, 1]
-end
-
+
 # Extracts a sub-data set using indices
 function ExtractDataSets(V::AbstractGeneralGrid, iLon, iLat, iDepth)
 

test/test_data_import.jl

@@ -78,3 +78,10 @@ data_Image          =   Screenshot_To_UTMData(filename,Corner_LowerLeft, Corner_
 @test data_Image.EW.val[22] ≈ 0.42424242424242425
 @test data_Image.NS.val[22] ≈ 48.666666666666664
 @test Value(data_Image.depth[22]) ≈ -15.0m
+
+# test the import of xml files from ISC
+# the search criteria are set in a way that only one event should be found
+download_data("http://www.isc.ac.uk/cgi-bin/web-db-run?request=COLLECTED&req_agcy=ISC-EHB&out_format=QuakeML&ctr_lat=&ctr_lon=&radius=&max_dist_units=deg&searchshape=RECT&top_lat=49&bot_lat=37&left_lon=4&right_lon=20&srn=&grn=&start_year=2000&start_month=1&start_day=01&start_time=00%3A00%3A00&end_year=2005&end_month=12&end_day=31&end_time=00%3A00%3A00&min_dep=&max_dep=&min_mag=5.8&max_mag=&req_mag_type=Any&req_mag_agcy=Any&min_def=&max_def=&include_magnitudes=on&include_links=on&include_headers=on&include_comments=on&table_owner=iscehb","ISCTest.xml")
+Data_ISC = GetLonLatDepthMag_QuakeML("ISCTest.xml");
+@test Value(Data_ISC.depth[1])==-13.0km
+@test Data_ISC.fields.Magnitude[1]==5.8

test/test_utils.jl

@@ -10,6 +10,16 @@ Vx1,Vy1,Vz1     =   Data1*3,Data1*4,Data1*5
 Data_set2D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data1,LonData1=Lon, Velocity=(Vx1,Vy1,Vz1)))  
 @test_throws ErrorException CrossSection(Data_set2D, Depth_level=-10)
 
+# Test interpolation of depth to a given cartesian XY-plane
+x = 11:19
+y = 31:39
+plane1 = InterpolateDataFields2D(Data_set2D, x, y)
+proj   = ProjectionPoint()
+plane2 = InterpolateDataFields2D(Data_set2D, proj, x, y)
+
+@test plane1 == plane2
+@test all(==(-50e0), plane1)
+
 # Create 3D volume with some fake data
 Lon,Lat,Depth   =   LonLatDepthGrid(10:20,30:40,(-300:25:0)km);
 Data            =   Depth*2;                # some data