RyallBFPcode.nb

Code by Vince Cronin

This code is adapted from the code presented in the following resource:
Paláncz, B., Lovas, T., Molnár, B., 2013, Plane fitting to point cloud via Gröbner basis:  Wolfram Library Archive, accessed 2020 Nov. 2 via https://library.wolfram.com/infocenter/ID/8769/

Version of December 19, 2020

Clear[locXYZ,ℕ]

unitVect3D[vect_]:={(vect[[1]]/Norm[vect]),(vect[[2]]/Norm[vect]),(vect[[3]]/Norm[vect])};

Data are from Table 2 of
Ryall, A., VanWormer, J.D., and Jones, A.E., 1968, Triggering of microearthquakes by Earth tides and other features of the Truckee, California, earthquake sequence of September 1966:  Bulletin of the Seismological Society of America, v. 58, no. 1, p. 215-248.

inData=Import["/Users/vince/Desktop/Dec19-BFP-R/RyallNAD27UTM.xlsx"];

eqData=Flatten[inData,1];

noRecords=Length[eqData]

a = Max[Table[eqData[[i, 2]], {i, 1, noRecords}]]

b = Min[Table[eqData[[i, 2]], {i, 1, noRecords}]]

d = Max[Table[eqData[[i, 3]], {i, 1, noRecords}]]

e = Min[Table[eqData[[i, 3]], {i, 1, noRecords}]]

The UTM coordinates of the lower - left corner (below) are somewhat arbitrary, but were chosen to coincide with the lower-left corner of “panel 52” of the hillshade mosaic made from LiDAR data from the 2014 National Forest dataset accessed on OpenTopography.org.  That is the southwest corner of the mosaic.

xllcorner= 739043;

yllcorner= 4358540;

zllcorner=0;

locXYZ=Table[{(Round[eqData[[i,2]]]-xllcorner),(Round[eqData[[i,3]]]-yllcorner),(zllcorner-(1000*eqData[[i,4]]))},{i,1,noRecords}];

ℕ=Length[locXYZ]

a = Max[Table[locXYZ[[i, 1]], {i, 1, noRecords}]]

b = Min[Table[locXYZ[[i, 1]], {i, 1, noRecords}]]

d = Max[Table[locXYZ[[i, 2]], {i, 1, noRecords}]]

e = Min[Table[locXYZ[[i, 2]], {i, 1, noRecords}]]

f = Max[Table[locXYZ[[i, 3]], {i, 1, noRecords}]]

g = Min[Table[locXYZ[[i, 3]], {i, 1, noRecords}]]

p8=ListPointPlot3D[locXYZ,PlotStyle->Blue];

Show[p8,ViewPoint->Above]

Show[p8, ViewPoint -> {0, 0, ∞}]

Solution via SVD

c0=Mean[locXYZ];

N[%]

Cc0=Map[#-c0&,locXYZ];

AbsoluteTiming[{u,w,v}=SingularValueDecomposition[Cc0];];

ν=Last[Transpose[v]];

({𝓍,𝓎,𝓏}-c0).ν==0//Simplify;

solnSVD=NSolve[%,𝓏]//Expand;

N[%]

p4=Plot3D[𝓏/.%[[1]],{𝓍,b,a},{𝓎,d,e},PlotRange->{g,f},BoxRatios->{1,1,1}, ClippingStyle->None,Mesh->None,ColorFunction->"SandyTerrain", PlotStyle->Opacity[0.3]];

Show[{p8,p4},ViewPoint->Above]

Show[{p8,p4}, ViewPoint -> {0, 0, ∞}]

Solution via PCA

AbsoluteTiming[A=Covariance[Cc0];MatrixForm[A]];

{nx,ny,nz}=Eigenvectors[A][[3]];

sold=Solve [{nx,ny,nz}.Total[locXYZ]-𝒹 ℕ==0,𝒹]//Flatten;

solnPCA={-(nx/nz),-(ny/nz),(𝒹/nz)}/.sold;

N[%]

Doodling

α=solnPCA[[1]];

N[%]

β=solnPCA[[2]];

N[%]

γ=solnPCA[[3]];

N[%]

meanX=c0[[1]];

meanY=c0[[2]];

meanZ=c0[[3]];

Tests

meanX;

x=(γ+(β*meanY)-meanZ)/(-α);

meanY;

y=(γ+(α*meanX)-meanZ)/(-β);

meanZ;

z=γ+(α*meanX)+(β*meanY);

What is the RHR strike and dip angle of the plane?

a2=Max[Table[locXYZ[[i,1]],{i,1,noRecords}]];

b2=Min[Table[locXYZ[[i,1]],{i,1,noRecords}]];

midX=(a2+b2)/2;

d2=Max[Table[locXYZ[[i,2]],{i,1,noRecords}]];

e2=Min[Table[locXYZ[[i,2]],{i,1,noRecords}]];

midY=(d2+e2)/2;

trialX=(a2+b2)/2;

trialY=(γ+(α*trialX)-meanZ)/(-β);

lowZ=meanZ-3;

lowY=(γ+(α*trialX)-lowZ)/(-β);

onStrike={trialX,trialY,meanZ};

meanVector={meanX,meanY,meanZ};

belowStrike={trialX,lowY,lowZ};

strikeVect1=onStrike-meanVector;

obliqueVect=belowStrike-meanVector;

normal1=unitVect3D[Cross[strikeVect1,obliqueVect]];

normal2=-normal1;

refNormal=If[(normal1[[3]]>normal2[[3]]),normal1,If[(normal1[[3]]<normal2[[3]]),normal2,normal1]];

up={0,0,1};

dipAngle=VectorAngle[refNormal,up]*(180/π);

projectedVector={refNormal[[1]],refNormal[[2]],0};

unitProjVector=unitVect3D[projectedVector];

north={0,1,0};

fromNorth=VectorAngle[north,unitProjVector]*(180/π);

dipAzimuth=If[(unitProjVector[[1]]≥0),fromNorth,(360-fromNorth)];

N[%]

rhrGridStrike=If[(dipAzimuth≥90),(dipAzimuth-90),(dipAzimuth+270)];

N[%]

Adjustment for input geographic coordinates based on true north to the relevant UTM grid north

zoneMeridian=-123;

nwLat=39.5801; nwLong=-120.2732;
swLat=39.3375; swLong=-120.2727;
seLat=39.3375; seLong=-119.9875;
neLat=39.5801;neLong=-119.9878;

midLat=((nwLat-swLat)/2)+swLat;

N[%]

midLongN=((neLong-nwLong)/2)+nwLong;

midLongS=((seLong-swLong)/2)+swLong;

midLong=((midLongS-midLongN)/2)+midLongN;

N[%]

The gridConvergence is the angle between true north and grid north at/near the center of the map area.  If grid north is found by a clockwise rotation from true north, the sign of the gridConvergence is positive;  otherwise, it is negative.  
Imagine you have a horizontal vector pointing toward true north (azimuth 0°).  If the longitude of the point-of-action (or the middle of the map area) is greater than the longitude of the central meridian of the UTM zone (i.e., > -123 for Zone 10_), the azimuth of that vector in the UTM grid at that point-of-action will be the true azimuth minus the gridNorthAdjustment.

gridConvergence=ArcTan[Tan[(midLong-zoneMeridian)Degree]*Sin[midLat Degree]]*(180/π);

N[%]

rhrStrike=(rhrGridStrike+gridConvergence);

N[%]

N[dipAngle]

What is the location (UTM easting, northing, depth) of a point along the best - fit plane?

recomputedVector={(meanX+xllcorner),(meanY+yllcorner),meanZ};

meanUTMeasting=recomputedVector[[1]];

N[%]

meanUTMnorthing=AccountingForm[recomputedVector[[2]]];

N[%,7]

meanDepth=recomputedVector[[3]];

N[%]

ClearAll[meanDepth,meanUTMnorthing,meanUTMeasting,recomputedVector,rhrStrike,gridConvergence,rhrGridStrike,dipAzimuth,fromNorth,unitProjVector,refNormal,normal2,normal1,obliqueVect,eqData,noRecords,p8,c0,Cc0,ν,p4]

Summary — V4

The best - fit plane through the composite list of aftershock focal locations has a right - hand - rule strike of 206.123° and a dip angle of 89.6472° toward the NW.  In NAD27, that plane passes through a point at latitude 39.4101218708°, longitude -120.1740925339°, depth 4.63718 km subsea (UTM Zone 10N 743307 mE 4365892 mN).

Geographic coordinates in NAD27 computed using the NGS Coordinate Conversion and Transformation Tool at https://www.ngs.noaa.gov/NCAT/index.xhtml