From 72a61ca06d97e9f8afac661b3ba3c714b4a9d812 Mon Sep 17 00:00:00 2001 From: Albert de Montserrat Date: Thu, 29 Feb 2024 12:26:44 +0100 Subject: [PATCH 1/3] 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 2/3] 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 3/3] 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)