Hello,
Is there a way to contour or use Sgrayplot to create a density plot of this Kaverina data plot? Meshgrid seemed to fail. The scatter plot is fine.
Thanks
Lester
Hello,
What is the data ? For arbitrary domain you should have a triangular mesh of the domain, then Fec and contour2dm can probably do what you want.
S.
The data are X, Y, and Z=depth. X and Y are coordinates in the Kaverina plot (kav_x, Kav_y).The data do not extend beyond the base plot.
The Kaverina plot is a way of plotting earthquake data and getting a visualisation of the type of mechanism (rupture type) from the seismic moment tensor.
How would a triangular mesh work for this type of plot? I gridded the data in Surfer and areas beyond the plot region are replaced with NaNs. If it is helpful, I can e-mail the CSV file?
The inputs all have a size of 3606x1.
Lester
For a more classical ternary plot, which amounts to consider a function f(x,y) on the domain
\{x\geq 0, y\geq 0,x+y\leq 1\}
and represent either a surface or a pseudo-color map, I have written the ternmesh function below, which yields the triangular mesh of such a triangle, which is easily transformed to the classical equilateral triangle of ternary plots. Then f(x,y) is evaluated on original coordinates and the actual abscissae of the triangles of the plot are obtained by a matrix product (see below).
Below I have considered a use case where f(x,y) is a density of random points smoothed by a Gaussian kernel. I think that you can adapt the code to your needs:
function [x,y,tri] = ternmesh(n)
[i,j] = meshgrid(1:n);
[x,y] = meshgrid(linspace(0,1,n));
k = (i-1)*n+j;
flag = (i+j)<=n+1;
x = x(flag);
y = y(flag);
k(flag)=1:n*(n+1)/2;
k(~flag)=%nan;
tri = [k(1:$-1,1:$-1)(:) k(2:$,1:$-1)(:) k(1:$-1,2:$)(:)
k(2:$,2:$)(:) k(1:$-1,2:$)(:) k(2:$,1:$-1)(:)];
tri(or(isnan(tri),2),:) = [];
endfunction
[x,y,tri]=ternmesh(50);
id = 1:size(tri,1);
z=0;
n=0;
for i=1:20
a=rand();
b=rand();
v=rand();
if a+b <= 1
n = n+1;
sigma = 0.1;
z = z + exp(-((x-a).^2+(y-b).^2)/sigma^2)/(sqrt(2*%pi)*sigma)
end
end
z = z/n;
colormap(parula(128))
xy = [1 1/2;0 sqrt(3)/2]*[x(:), y(:)]';
fec(xy(1,:),xy(2,:),[id' tri id'],z)
isoview("on")
colorbar
Hi Stephane,
The Kaverina diagram is similar to a ternary plot, but is a projected coordinate system more like a “segment” of the Lambert Equal Area projection. I guess that gridding the data should be over a square/rectangular domain that contains the data itself (as Surfer does).
The Kaverina digram (for X and Y) is defined by the following function:
function [x, y] = kave(plungt, plungb, plungp)
// Computes x and y for the Kaverina diagram
zt = sin(deg2rad(plungt));
zb = sin(deg2rad(plungb));
zp = sin(deg2rad(plungp));
L = 2 * sin(0.5 * acos((zt + zb + zp) / sqrt(3)));
N = sqrt(2 * ((zb - zp)^2 + (zb - zt)^2 + (zt - zp)^2));
x = sqrt(3) * (L / N) * (zt - zp);
y = (L / N) * (2 * zb - zp - zt);
endfunction
The Kaverina digram plotting:
function plot_P_axis()
X = zeros(1,101);
Y = zeros(1,101);
for a = 0:100
P = asin(sqrt(a/100.0)) * (180/%pi);
B = 0.0;
T = asin(sqrt(1-(a/100.0))) * (180/%pi);
[X(a+1), Y(a+1)] = kave(T, B, P);
end
tickx = zeros(1,10);
ticky = zeros(1,10);
for i = 0:9
[tickx(i+1), ticky(i+1)] = kave(90 - 10*i, 0, 10*i);
end
plot(X, Y, 'r-', 'LineWidth', 2); // P-axis
//scatter(tickx, ticky); // P-points
scatter(tickx, ticky-0.018, 'marker', 20); // P-points
gce.children.mark_size=10
gce.children.mark_foreground=-1
gce.children.thickness=2
for i = 0:9
xstring(tickx(i+1)-0.023, ticky(i+1)-0.1, string(10*i));
end
xstring(-0.15, -0.76, 'P axis plunge');
endfunction
function plot_B_axis()
X = zeros(1,101);
Y = zeros(1,101);
for a = 0:100
B = asin(sqrt(a/100.0)) * (180/%pi);
P = 0.0;
T = asin(sqrt(1-(a/100.0))) * (180/%pi);
[X(a+1), Y(a+1)] = kave(T, B, P);
end
tickx = zeros(1,10);
ticky = zeros(1,10);
for i = 0:9
[tickx(i+1), ticky(i+1)] = kave(0, 10*i, 90 - 10*i);
end
plot(X, Y, 'r-', 'LineWidth', 2); //T-axis
//scatter(tickx, ticky); // B-points
scatter(tickx-0.018, ticky, 'marker', 19); // B-points
gce.children.mark_size=10
gce.children.mark_foreground=-1
gce.children.thickness=2
for i = 0:9
xstring(tickx(i+1)-0.04, ticky(i+1), string(10*i));
end
xstring(-0.7, 0.2, 'B axis plunge');
gce.font_angle=-65
endfunction
function plot_T_axis()
X = zeros(1,101);
Y = zeros(1,101);
for a = 0:100
T = asin(sqrt(a/100.0)) * (180/%pi);
B = 0.0;
P = asin(sqrt(1-(a/100.0))) * (180/%pi);
[X(a+1), Y(a+1)] = kave(T, B, P);
end
tickx = zeros(1,10);
ticky = zeros(1,10);
for i = 0:9
[tickx(i+1), ticky(i+1)] = kave(10*i, 90 - 10*i, 0);
end
//plot(X, Y, 'k-', 'LineWidth', 2); // P-axis
scatter(tickx+0.018, ticky, 'marker', 19);// T-points
gce.children.mark_size=10
gce.children.mark_foreground=-1
gce.children.thickness=2
//scatter(tickx + 0.025, ticky);// T-points (original offset from axis)
for i = 0:9
xstring(tickx(i+1)+0.01, ticky(i+1)-0.01,string(10*i));
end
xstring(0.61, 0.45, 'T axis plunge');
gce.font_angle=65
X = zeros(1,101);
Y = zeros(1,101);
for a = 0:100
P = asin(sqrt(a/100.0)) * (180/%pi);
T = 0.0;
B = asin(sqrt(1-(a/100.0))) * (180/%pi);
[X(a+1), Y(a+1)] = kave(T, B, P);
end
plot(X, Y, 'r-', 'LineWidth', 2); // B-axis
endfunction
Plot the Kaverina base map
plot_P_axis()
plot_B_axis()
plot_T_axis()
xstring(-1.0, -0.5, 'Normal');
xstring(0.9, -0.5, 'Reverse');
xstring(-0.08, 1, 'Strike-Slip');
f=gcf()
f.figure_size=[1024, 1024]
a = gca();
a.isoview='on'
a.box='off'
a.axes_visible=['off', 'off', 'off']
a.data_bounds = [-1.2,-1.2;1.2,1.2]
Lester
Reading the Kaverina paper (https://academic.oup.com/gji/article-pdf/125/1/249/5891249/125-1-249.pdf) seems to be a must to understand (I deleted my previous answer which contained questions that are answered in the paper). I also guess that you started by converting a Python code ?
Anyway, can you post the part of your code that plots the actual data on the diagram ?
S.
Hi Stephane,
Yes I converted the Python code, and tried adding more comments to explain what is going on.
Running the Kaverina baseplot function and then a simple scatter plot gives the plot. Normally one would define a grid domain slightly larger than the data region and make a grid like that such that areas with no data are NaN.
I have sent you the data files to check via mail. That way you can see the logic and make sure it is working correctly. The Python code has not changed that much.
Lester
Hello,
After vectorizing your kave() function:
function [x, y] = kave(plungt, plungb, plungp)
// Computes x and y for the Kaverina diagram
zt = sin(deg2rad(plungt));
zb = sin(deg2rad(plungb));
zp = sin(deg2rad(plungp));
L = 2 * sin(0.5 * acos((zt + zb + zp) / sqrt(3)));
N = sqrt(2 * ((zb - zp).^2 + (zb - zt).^2 + (zt - zp).^2));
x = sqrt(3) * (L ./ N) .* (zt - zp);
y = (L ./ N) .* (2 * zb - zp - zt);
endfunction
and careful reading of the Kaverina paper, I have concluded that a trivial parametrization of the (T,P,B) domain is
\begin{align}
T&=\arcsin(\cos \theta\cos \varphi),\\
B&=\arcsin(\sin\theta\cos \varphi),\\
P&=\arcsin(\sin\varphi),
\end{align}
for (\theta,\phi) in S=[0,\pi/2]\times[0,\pi/2]. The square S cannot be meshed uniformly because all points having \phi=\pi/2 gather for all \theta on the north pole at T=0,B=0,P=\pi/2. That’s where my trimesh() function is useful as it can be used to mesh S with one node for the pole and a linear progression of the number of nodes for each \varphi when it goes from \pi/2 to 0, by using a simple non-linear transformation:
[u,v,tri] = ternmesh(30);
th = u./(1-v)*%pi/2;
th(v==1) = 0;
ph = v*%pi/2;
[x,y] = kave(asin(cos(th).*cos(ph))*180/%pi,...
asin(sin(th).*cos(ph))*180/%pi,...
ph*180/%pi);
it yields a quasi-uniform mesh in terms of projected (x,y):
scatter(x,y)
Then, it is quite easy to represent any scalar field, for example the kernel density of random points:
z = 0;
n = 100;
sigma = 0.1;
for i=1:n
th = rand()*%pi/2;
ph = rand()*%pi/2;
[a,b] = kave(asin(cos(th).*cos(ph))*180/%pi,...
asin(sin(th).*cos(ph))*180/%pi,...
ph*180/%pi);
z = z + exp(-((x-a).^2+(y-b).^2)/sigma^2)/(sqrt(2*%pi)*sigma)/n
end
id = 1:size(tri,1);
colormap(parula(128))
fec(x,y,[id' tri id'],z)
isoview("on")
colorbar
I think that’s what you needed, and at the end it was just a matter of correctly meshing a eigth of a sphere!
S.
Many thanks for the guidance and solution Stephane. You solved the problem!
Lester
1 Like