From f6893499e397c6b0604e5a214e0ab6d1bc9165a9 Mon Sep 17 00:00:00 2001 From: mthielma Date: Wed, 28 Feb 2024 19:56:06 +0100 Subject: [PATCH 1/5] added a function to import lon/lat/depth/mag from QuakeML files from ISC --- Project.toml | 3 ++- src/data_import.jl | 42 ++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 42 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 2334aa46..c35ce00d 100644 --- a/Project.toml +++ b/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"] diff --git a/src/data_import.jl b/src/data_import.jl index 247ac27f..f29d1a60 100644 --- a/src/data_import.jl +++ b/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,depth,(Magnitude=mag,Depth=depth)); + return Data_ISC end From 72a61ca06d97e9f8afac661b3ba3c714b4a9d812 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Thu, 29 Feb 2024 12:26:44 +0100 Subject: [PATCH 2/5] export `InterpolateTopographyOnPlane` as `InterpolateDataFields2D` --- src/utils.jl | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 4f8dc2cc..7f7ae82e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -2,7 +2,7 @@ export meshgrid, CrossSection, CrossSectionVolume, CrossSectionSurface, CrossSectionPoints, ExtractSubvolume, SubtractHorizontalMean export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap -export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane +export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields export RotateTranslateScale export DrapeOnTopo, LithostaticPressure! export FlattenCrossSection @@ -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 InterpolateTopographyOnPlane(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) @@ -1327,18 +1344,6 @@ 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) From 9acf29d134d0747985636a27f4a5d9745802f14e Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Thu, 29 Feb 2024 12:27:01 +0100 Subject: [PATCH 3/5] add unit test for new `InterpolateDataFields2D` method --- test/test_utils.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_utils.jl b/test/test_utils.jl index 573fe49d..e2f28e83 100644 --- a/test/test_utils.jl +++ b/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 From 78c7be027c8b5ee1d434b1339553aef7f0088197 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> Date: Thu, 29 Feb 2024 13:50:54 +0100 Subject: [PATCH 4/5] fix function call --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index ecc22983..e45c4e71 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1182,7 +1182,7 @@ 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 InterpolateTopographyOnPlane(V::GeoData, proj, x, y) + return InterpolateDataFields2D(V::GeoData, proj, x, y) end function InterpolateDataFields2D(LonLat::GeoData, proj::ProjectionPoint, x::AbstractRange, y::AbstractRange) From c562f0e5ab2e251c847bcf3ff6b6827286065b8c Mon Sep 17 00:00:00 2001 From: mthielma Date: Thu, 29 Feb 2024 14:08:32 +0100 Subject: [PATCH 5/5] updated data_import.jl - added conversion to km when creating the GeoData object - added test --- src/data_import.jl | 2 +- test/test_data_import.jl | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/data_import.jl b/src/data_import.jl index f29d1a60..16576256 100644 --- a/src/data_import.jl +++ b/src/data_import.jl @@ -348,6 +348,6 @@ function GetLonLatDepthMag_QuakeML(filename::String) 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,depth,(Magnitude=mag,Depth=depth)); + Data_ISC = GeoData(lon,lat,-1*depth/1e3,(Magnitude=mag,Depth=-1*depth/1e3*km)); return Data_ISC end diff --git a/test/test_data_import.jl b/test/test_data_import.jl index eafc94f0..883c7152 100644 --- a/test/test_data_import.jl +++ b/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 \ No newline at end of file