gps-strain-calculator-mathematica.nb

Code to determine the 2-D (horizontal) instantaneous or infinitesimal strain from velocity or displacement data from 3 adjacent-and-non-colinear GNSS stations

Coded in Mathematica 8 by Vince Cronin with help from Phil Resor, circa 2012.  This version is coded in Mathematica 11, and was revised 13 October 2018.
The most current version of this code is accessible via 
http://croninprojects.org/Vince/Codes/gps_strain_calculator_mathematica.nb

Input

Data for this analysis can be obtained from many sources, but you must be certain that the velocities for all three sites are expressed with respect to the same reference frame (e.g., NAM08 or IGS08).
The Plate Boundary Observatory GPS data are available online through the good work of UNAVCO via
http://www.unavco.org/instrumentation/networks/status/pbo/gps
Each row of the input data contains the following data types :  GPS station number, longitude (west is negative), latitude (south is negative), E-W velocity (mm/yr), N-S velocity (mm/yr), E-W velocity uncertainty (mm/yr), N-S velocity uncertainty (mm/yr).  
The following is an example of a valid input matrix.

inData=(({{P574, -117.63389, 34.28677, -19.23, 19.61, 0.16, 0.19}, {P577, -117.31891, 34.30461, -14.51, 16.15, 0.9, 0.09}, {EWPP, -117.52559, 34.10420, -21.83, 22.56, 0.06, 0.07}}));

Replace the data in the inData matrix below with your values.

inData=(({{P574, -117.63389, 34.28677, -19.23, 19.61, 0.16, 0.19}, {EWPP, -117.52559, 34.10420, -21.83, 22.56, 0.06, 0.07}, {P577, -117.31891, 34.30461, -14.51, 16.15, 0.9, 0.09}}));

UTM Central Meridian

Input the central meridian that is spatially appropriate for the three PBO sites used in this analysis.  If the three points involve two adjacent UTM zones, input the central meridian of the zone to the west.

({{Longitude Range, Central Meridian}, {-138 to -132, -135}, {-132 to -126, -129}, {-126 to -120, -123}, {-120 to -114, -117}, {-114 to -108, -111}, {-108 to -102, -105}, {-102 to -96, -99}})

Replace the value for longCentMeridian below with your value.

longCentMeridian=-117;

The results are in the output section at the end of this notebook

User - Defined Functions

The purpose of the following function is to convert input data from geographic coordinates (latitude and longitude) to UTM using the UTM84 datum.  Online converters are available via https://www.ngs.noaa.gov/NCAT/,  http://www.rcn.montana.edu/Resources/Converter.aspx, 
https://awsm-tools.com/geo/utm-to-geographic, and elsewhere.  The output data are

geog2UTM[inputLat_,inputLong_]:=Module[{c42,c43,c44,c45,c46,c47,c48,c49,c50,c51,c53,c54,c56,c57,c58,c59,c60,c61,c62,c63,c65,c66,c67,c68,c69,c71,c72,c73,c74,c75,c78,c79,c80,c81,outData},c42=inputLat (π/180);c43=inputLong (π/180);c44=(inputLong+180)/6;c45=If[(c44-IntegerPart[c44])>0,(IntegerPart[c44]+1),IntegerPart[c44]];(* c45 is the central meridian used to compute the location *)
c46=-183+(6*c45);c47=c46 (π/180);c48=If[(c46==-177),177,c46-6];c49=c48 (π/180);c50=(c48+180)/6;c51=If[(c50-IntegerPart[c50])>0,IntegerPart[c50]+1,IntegerPart[c50]];c53=6378137;c54=6356752.3142;c56=0.9996;c57=(Sqrt[(1-(c54^2)/(c53^2))]);c58=((((c57*c53)/c54))^2);c59=(c53-c54)/(c53+c54);c60=(c53*(1-(c57^2)))/((((1-((c57^2)*((Sin[c42])^2))))^1.5));c61=c53/(Sqrt[1-((c57^2)*((Sin[c42])^2))]);c62=c43-c47;c63=If[(c48==177),(Abs[c43]-c49),(c43-c49)];c65=c42*(1-((c57^2)/4)-((3*(c57^4))/64)-((5*(c57^6))/256));c66=Sin[2*c42]*(((3*(c57^2))/8)+((3*(c57^4))/32)+((45*(c57^6))/1024));c67=Sin[4*c42]*(((15*(c57^4))/256)+((45*(c57^6))/1024));c68=Sin[6*c42]*((35*(c57^6))/3072);c69=c53*(c65-c66+c67-c68);c71=c69*c56;c72=c56*c61*Sin[2*c42]/4;c73=(c56*c61*Sin[c42]*(((Cos[c42])^3))/24)*(5-(((Tan[c42])^2))+(9*c58*(((Cos[c42])^2)))+(4*(c58^2)*((Cos[c42])^4)));c74=c56*c61*Cos[c42];c75=(c56*c61*(((Cos[c42])^3))/6)*(1-((Tan[c42])^2)+(c58*((Cos[c42])^2)));c78=c71+(c72*(c62^2))+(c73*(c62^4));(* c78 is the directly computed northing*)
c79=500000+(c74*c62)+(c75*(c62^3));(* c79 is the directly computed easting*)
c80=c71+(c72*(c63^2))+(c73*(c63^4));(* c80 is the pseudo-northing based on the central meridian to the west *)
c81=500000+(c74*c63)+(c75*(c63^3));(* c81 is the pseudo-easting based on the central meridian to the west *)
outData={c78,c79,c80,c81,c45};outData];

The purpose of the following function is to jointly consider the UTM coordinates of three sites that might be in two adjacent UTM zones and re - express their coordinates in terms of the west - most zone of the set.

triGroupUTM[site1_,site2_,site3_]:=Module[{d83,utmZones,e83,c83,c87,d87,e87,c88,d88,e88,finalEastings,finalNorthings,answer},d83=(site1[[5]]+site2[[5]]+site3[[5]])/3;utmZones={site1[[5]],site2[[5]],site3[[5]]};e83=StandardDeviation[utmZones];c83=If[e83>5,60,Floor[d83]];c87=If[site1[[5]]==c83,site1[[2]],site1[[4]]];d87=If[site2[[5]]==c83,site2[[2]],site2[[4]]];e87=If[site3[[5]]==c83,site3[[2]],site3[[4]]];
c88=If[site1[[5]]==c83,site1[[1]],site1[[3]]];d88=If[site2[[5]]==c83,site2[[1]],site2[[3]]];e88=If[site3[[5]]==c83,site3[[1]],site3[[3]]];finalEastings={c87,d87,e87};finalNorthings={c88,d88,e88};answer={finalEastings,finalNorthings};answer
];

The purpose of the following function is to create a plot of the infinitesimal horizontal strain ellipse, exaggerated by a factor of one million.

StrainEllipse[maxHextension_,minHextension_,azimuthS1H_,azimuthS2H_]:=Module[{a,circleData,circle,a2,b,line1,line2,redLine,blueLine,ellipseData,zMatrix,rotatedEllipseData,ellipsePlot},a=1;circleData=Table[{a Cos[i],a Sin[i]},{i,0,2π,(π/50)}];circle=Graphics[Line[circleData],PlotRange->2,Axes->True,LabelStyle->{FontFamily->"Arial",10,GrayLevel[0]},AxesLabel->{HoldForm["East"],HoldForm["North"]},PlotLabel->HoldForm["Strain ellipse ((10^6) exaggeration)"],ImageSize->Large];a2=1+(maxHextension*0.001);b=1+(minHextension*0.001);
line1=Line[{{a2*Sin[azimuthS1H Degree],a2*Cos[azimuthS1H Degree]},{-a2*Sin[azimuthS1H Degree],-a2*Cos[azimuthS1H Degree]}}];line2=Line[{{b*Sin[azimuthS2H Degree],b*Cos[azimuthS2H Degree]},{-b*Sin[azimuthS2H Degree],-b*Cos[azimuthS2H Degree]}}];
redLine=Graphics[{Red,line1}];
blueLine=Graphics[{Blue,line2}];
ellipseData=Table[{b Cos[i],a2 Sin[i]},{i,0,2π,(π/50)}];
zMatrix=(({{Cos[-azimuthS1H Degree], -Sin[-azimuthS1H Degree]}, {Sin[-azimuthS1H Degree], Cos[-azimuthS1H Degree]}}));
rotatedEllipseData=Flatten[Table[{zMatrix.ellipseData[[i]]},{i,1,101}],1];ellipsePlot=Graphics[Line[rotatedEllipseData],PlotRange->2,Axes->True,LabelStyle->{FontFamily->"Arial",10,GrayLevel[0]},AxesLabel->{HoldForm["East"],HoldForm["North"]},PlotLabel->HoldForm["Strain ellipse ((10^6) exaggeration)"],ImageSize->Large];Show[circle,redLine,blueLine,ellipsePlot]];

The purpose of the following function is to determine the azimuth and speed (magnitude) of a total velocity vector, given the east and north velocity components.  The units output are degrees of azimuth and total speed in the same units as the input velocity values (e.g., mm/yr).

azimuthNSpeed[eastVelocity_,northVelocity_]:=Module[{unitNorthVector,totalVelVector,totalSpeed,unitTotalVelVector,angle,azTotalVelVector,answer},unitNorthVector={0,1};totalVelVector={eastVelocity,northVelocity};totalSpeed=Norm[totalVelVector];unitTotalVelVector={eastVelocity/totalSpeed,northVelocity/totalSpeed};angle=ArcCos[unitTotalVelVector.unitNorthVector]((180/π));azTotalVelVector=If[(unitTotalVelVector[[1]]<0),360-angle,angle];answer={azTotalVelVector,totalSpeed};answer];

Calculations

The steps referenced below are described in the document Algorithm for computing infinitesimal strain rate between three non-colinear GPS stations, given their N-S and E-W velocities, with a worked example (Vince Cronin and Phil Resor, 2012; revised 2018).  The most current version of that document is available online via http://croninprojects.org/Vince/Geodesy/Current-Strain-Calculators/TriangleStrainAlgorithm.pdf

Step 1:  Convert latitude & longitude to UTM

site1=inData[[1,1]];site2=inData[[2,1]];site3=inData[[3,1]];

site1UTM=geog2UTM[inData[[1,3]],inData[[1,2]]];

site2UTM=geog2UTM[inData[[2,3]],inData[[2,2]]];

site3UTM=geog2UTM[inData[[3,3]],inData[[3,2]]];

finalUTMcoord=triGroupUTM[site1UTM,site2UTM,site3UTM];

eastings={finalUTMcoord[[1,1]],finalUTMcoord[[1,2]],finalUTMcoord[[1,3]]};

northings={finalUTMcoord[[2,1]],finalUTMcoord[[2,2]],finalUTMcoord[[2,3]]};

Step 2:  Find the center of the triangle defined by the three sites

meanX=(1/3)(eastings[[1]]+eastings[[2]]+eastings[[3]]);

meanY=(1/3)(northings[[1]]+northings[[2]]+northings[[3]]);

meanLat=(inData[[1,3]]+inData[[2,3]]+inData[[3,3]])/3;

meanLong=(inData[[1,2]]+inData[[2,2]]+inData[[3,2]])/3;

Steps 3 & 4:  Define the initial locations relative to the center of the triangle

m1=(({{eastings[[1]]-meanX, northings[[1]]-meanY}, {eastings[[2]]-meanX, northings[[2]]-meanY}, {eastings[[3]]-meanX, northings[[3]]-meanY}}));

The first row of matrix m1 (above) contains the east - west and north - south coordinates of site 1 in a coordinate system whose origin is the center of the triangle.  Units are meters, and the coordinate grid is aligned with the local UTM grid.

Step 5

m2=(({{1, 0, -m1[[1,2]], m1[[1,1]], m1[[1,2]], 0}, {0, 1, m1[[1,1]], 0, m1[[1,1]], m1[[1,2]]}, {1, 0, -m1[[2,2]], m1[[2,1]], m1[[2,2]], 0}, {0, 1, m1[[2,1]], 0, m1[[2,1]], m1[[2,2]]}, {1, 0, -m1[[3,2]], m1[[3,1]], m1[[3,2]], 0}, {0, 1, m1[[3,1]], 0, m1[[3,1]], m1[[3,2]]}}));

Step 6:  Invert matrix m2 to form matrix m3

m3=Inverse[m2];

Step 7:  Create a column matrix called m4 that contains the velocity data for the three GPS sites (m/yr)

The input data provides the magnitude of east-west and north-south components of the site velocity vector in a geographic coordinate system.  When these site data are converted to UTM coordinates, the divergence between the grid north and true north can be significant, so the site velocity vectors must be transformed into the UTM coordinate system using the grid convergence angle at each site.

Angular difference between grid north and true north for each site

Answers are given in degrees.  Positive values indicate that true north is a positive (anti-clockwise) rotation from grid north.  Negative values indicate that true north is a negative (clockwise) rotation from grid north.

site1GridConvergence=ArcTan[Tan[(inData[[1,2]]-longCentMeridian) Degree] Sin[inData[[1,3]] Degree]](180/π);

site2GridConvergence=ArcTan[Tan[(inData[[2,2]]-longCentMeridian) Degree] Sin[inData[[2,3]] Degree]](180/π);

site3GridConvergence=ArcTan[Tan[(inData[[3,2]]-longCentMeridian) Degree] Sin[inData[[3,3]] Degree]](180/π);

The grid convergence at the centroid of the triangle, represented by centroidGridConvergence, is given in degrees.  Positive values indicate that true north is a positive (anti-clockwise) rotation from grid north.  Negative values indicate that true north is a negative (clockwise) rotation from grid north.

centroidGridConvergence=ArcTan[Tan[(meanLong-longCentMeridian) Degree] Sin[meanLat Degree]](180/π);

transformSite1Vector=(({{Cos[site1GridConvergence Degree], Sin[site1GridConvergence Degree]}, {-Sin[site1GridConvergence Degree], Cos[site1GridConvergence Degree]}})).(({{inData[[1,4]]}, {inData[[1,5]]}}));

transformSite2Vector=(({{Cos[site2GridConvergence Degree], Sin[site2GridConvergence Degree]}, {-Sin[site2GridConvergence Degree], Cos[site2GridConvergence Degree]}})).(({{inData[[2,4]]}, {inData[[2,5]]}}));

transformSite3Vector=(({{Cos[site3GridConvergence Degree], Sin[site3GridConvergence Degree]}, {-Sin[site3GridConvergence Degree], Cos[site3GridConvergence Degree]}})).(({{inData[[3,4]]}, {inData[[3,5]]}}));

The division of each input value by 1000 achieves a conversion from mm/yr to m/yr (or from mm to m, if the problem is framed in terms of an average annual displacement)

m4=(({{transformSite1Vector[[1,1]]/1000}, {transformSite1Vector[[2,1]]/1000}, {transformSite2Vector[[1,1]]/1000}, {transformSite2Vector[[2,1]]/1000}, {transformSite3Vector[[1,1]]/1000}, {transformSite3Vector[[2,1]]/1000}}));

Step 8:  Multiply matrix m4 by matrix m3 to yield a 6x1 column matrix called m5

m5=m3.m4;

Step 9:  Find the infinitesimal translation vector for the triangle defined by the three GPS stations

translationVector={m5[[1,1]],m5[[2,1]]};

The coordinates of the translationVector (above) are {east-west velocity, north-south velocity} expressed in meters per year, relative to the stable North American reference frame (NAM08).  This can be considered either a velocity vector (m/yr) or an average annual displacement vector (m).

translationSpeed=(Sqrt[((m5[[1,1]])^2)+((m5[[2,1]])^2)]);

The translation speed represented by translationSpeed (above) is expressed in meters per year, relative to the stable North American reference frame (NAM08).  This can be considered either a speed (m/yr) or an average annual displacement (m).  Translation speed represented by translationSpeedMM (below) is expressed in millimeters/year.

translationSpeedMM=1000*translationSpeed;

unitTransVect={((translationVector[[1]])/(Norm[translationVector])),((translationVector[[2]])/(Norm[translationVector]))};

unitNorthVector={0,1};

angle=(ArcCos[unitTransVect.unitNorthVector]((180/π)))+centroidGridConvergence;

azTranslationVect=If[(unitTransVect[[1]]<0),360-angle,angle];

N[azTranslationVect];

The azimuth of the translation vector represented by azTranslationVect (above) is given in degrees, and is the horizontal angle measured clockwise from north to the direction the translation vector points.

Interlude:  Calculations related to creating the graphics associated with a typical GNSS-crustal-strain presentation,  assuming you already have a map graphic that shows the site locations

Note :  The site number (1, 2, or 3) is defined by the order in which the site data is listed in the input table (first row of input data is site 1)

Velocities in the North American reference frame (mm/yr)

The first number in each output pair is the azimuth in degrees;  the second number is the speed in millimeters per year (or average annual displacement in meters) relative to the NAM08 reference frame.

site1AzimuthNSpeed=azimuthNSpeed[inData[[1,4]],inData[[1,5]]];

site2AzimuthNSpeed=azimuthNSpeed[inData[[2,4]],inData[[2,5]]];

site3AzimuthNSpeed=azimuthNSpeed[inData[[3,4]],inData[[3,5]]];

The output of each of the preceding three lines of code are the site ID, azimuth, speed (mm/yr), east velocity, and north velocity relative to the NAM08 reference frame.

Velocities relative to the centroid of the GNSS station triangle (mm/yr), with the translation vector removed

site1NormAzNSpeed=azimuthNSpeed[(inData[[1,4]]-(1000*translationVector[[1]])),(inData[[1,5]]-(1000*translationVector[[2]]))];

site2NormAzNSpeed=azimuthNSpeed[(inData[[2,4]]-(1000*translationVector[[1]])),(inData[[2,5]]-(1000*translationVector[[2]]))];

site3NormAzNSpeed=azimuthNSpeed[(inData[[3,4]]-(1000*translationVector[[1]])),(inData[[3,5]]-(1000*translationVector[[2]]))];

The output of each of the preceding three lines of code are the site ID, azimuth, and speed (mm/yr) relative to the center of the triangle defined by the three PBO sites.

Distance between sites before displacement (km)

Distance between the first site and the second site (km)

sideLength12km=√((((m1[[1,1]]-m1[[2,1]]))^2)+(((m1[[1,2]]-m1[[2,2]]))^2))/1000;

Distance between the second site and the third site (km)

sideLength23km=√((((m1[[2,1]]-m1[[3,1]]))^2)+(((m1[[2,2]]-m1[[3,2]]))^2))/1000;

Distance between the first site and the third site (km)

sideLength13km=√((((m1[[1,1]]-m1[[3,1]]))^2)+(((m1[[1,2]]-m1[[3,2]]))^2))/1000;

Distance between sites after displacement (km)

Distance between the first site and the second site (km)

sideLengthEnd12km=√((((m10[[1,1]]-m10[[2,1]]))^2)+(((m10[[1,2]]-m10[[2,2]]))^2))/1000;

Distance between the second site and the third site (km)

sideLengthEnd23km=√((((m10[[2,1]]-m10[[3,1]]))^2)+(((m10[[2,2]]-m10[[3,2]]))^2))/1000;

Distance between the first site and the third site (km)

sideLengthEnd13km=√((((m10[[1,1]]-m10[[3,1]]))^2)+(((m10[[1,2]]-m10[[3,2]]))^2))/1000;

End Interlude

Step 10:  Find the value of the rotational velocity, Ω, in matrix m5: the value in (Subscript[m5, 31])

angularVelocityRad=m5[[3,1]];

angVelNanoRadYr=angularVelocityRad/(10^(-9));

Step 11:  Define the 2-D Lagrangian strain tensor ((Subscript[ε, ij])) and call it matrix m6

m6=(({{m5[[4,1]], m5[[5,1]]}, {m5[[5,1]], m5[[6,1]]}}));

Step 12:  Determine the orientations of the principal infinitesimal strain axes by finding the eigenvectors of the symmetrical Lagrangian infinitesimal strain tensor

eValues=Eigenvalues[m6];

e1=If[eValues[[1]]>eValues[[2]],eValues[[1]],eValues[[2]]];

e2=If[eValues[[1]]>eValues[[2]],eValues[[2]],eValues[[1]]];

e1-e2;

eVectors=Eigenvectors[m6];

s1V=If[eValues[[1]]>eValues[[2]],eVectors[[1]],eVectors[[2]]];

s1UV={((s1V[[1]])/(Norm[s1V])),((s1V[[2]])/(Norm[s1V]))};

s2V=If[eValues[[1]]>eValues[[2]],eVectors[[2]],eVectors[[1]]];

s2UV={((s2V[[1]])/(Norm[s2V])),((s2V[[2]])/(Norm[s2V]))};

angleS1HAxis=(ArcCos[unitNorthVector.s1UV]((180/π)))+centroidGridConvergence;

azimuthS1HAxis1=If[(s1UV[[1]]<0),(360-angleS1HAxis),angleS1HAxis];

azimuthS1HAxis2=If[azimuthS1HAxis1>180,azimuthS1HAxis1-180,azimuthS1HAxis1+180];

angleS2HAxis=(ArcCos[unitNorthVector.s2UV]((180/π)))+centroidGridConvergence;

azimuthS2HAxis1=If[s2UV[[1]]<0,360-angleS2HAxis,angleS2HAxis];

azimuthS2HAxis2=If[azimuthS2HAxis1>180,azimuthS2HAxis1-180,azimuthS2HAxis1+180];

Step 13:  Determine the magnitude of the maximum infinitesimal shear strain ((Subscript[γ, max])) at 45° to the maximum principal strain axis

maxShearStrain=2√(((((1/2)(m6[[1,1]]-m6[[2,2]])))^2)+(((m6[[1,2]]))^2));

Step 14:  Determine the areal strain between the GPS sites

circleArea=π;

s1=1+e1;

s2=1+e2;

ellipseArea=π×(s1)×(s2);

areaStrain=((ellipseArea-circleArea)/circleArea);

m6[[1,1]]+m6[[2,2]];

Step 15:  Identify/compute the second and third invariants of the strain tensor

firstInvariant=(m6[[1,1]]+m6[[2,2]]);

secondInvariant=((m6[[1,1]])(m6[[2,2]]))-(((m6[[1,2]]))^2);

thirdInvariant=Det[m6];

Step 16:  Compute the uncertainties

m7=(({{(1/((inData[[1,6]])^2)), 0, 0, 0, 0, 0}, {0, (1/((inData[[1,7]])^2)), 0, 0, 0, 0}, {0, 0, (1/((inData[[2,6]])^2)), 0, 0, 0}, {0, 0, 0, (1/((inData[[2,7]])^2)), 0, 0}, {0, 0, 0, 0, (1/((inData[[3,6]])^2)), 0}, {0, 0, 0, 0, 0, (1/((inData[[3,7]])^2))}}));

m8=Transpose[m2];

We use a Moore - Penrose matrix inverse (PseudoInverse) to compute M9, because the matrix resulting from m8.m7.m2 can be a poorly conditioned matrix.

m9=PseudoInverse[m8.m7.m2];

uncertTransE=(Sqrt[m9[[1,1]]]);

uncertTransN=(Sqrt[m9[[2,2]]]);

uncertRot=(Sqrt[m9[[3,3]]]);

uncertExx=(Sqrt[m9[[4,4]]]);

uncertExy=(Sqrt[m9[[5,5]]]);

uncertEyy=(Sqrt[m9[[6,6]]]);

summaryPlot=StrainEllipse[e1*(10^9),e2*(10^9),azimuthS1HAxis1,azimuthS2HAxis1];

Output

Distances between PBO sites before displacement

Distance between the first site and the second site (km)

{{site1,site2},sideLength12km}

Distance between the second site and the third site (km)

{{site2,site3},sideLength23km}

Distance between the first site and the third site (km)

{{site1,site3},sideLength13km}

Distances between PBO sites after displacement (km)

Distance between the first site and the second site (km)

{{site1,site2},sideLengthEnd12km}

Distance between the second site and the third site (km)

{{site2,site3},sideLengthEnd23km}

Distance between the first site and the third site (km)

{{site1,site3},sideLengthEnd13km}

Translation vector in the North American reference frame

Translation vector {east component, north component} where negative values are west and south, respectively.  Expressed in (mm/yr)

N[translationVector*1000]

Uncertainties in translation vector {east component, north component} expressed in (m/yr)

N[{uncertTransE,uncertTransN}]

Azimuth of translation vector (degrees)

N[azTranslationVect]

Translation speed (mm/yr)

N[translationSpeedMM]

Total velocities in the North American reference frame (mm/yr)

Site 1 azimuth and velocity in the North American reference frame (mm/yr)

{site1,(site1AzimuthNSpeed[[1]]+site1GridConvergence),site1AzimuthNSpeed[[2]]}

Site 2 azimuth and velocity in the North American reference frame (mm/yr)

{site2,(site2AzimuthNSpeed[[1]]+site2GridConvergence),site2AzimuthNSpeed[[2]]}

Site 3 azimuth and velocity in the North American reference frame (mm/yr)

{site3,(site3AzimuthNSpeed[[1]]+site3GridConvergence),site3AzimuthNSpeed[[2]]}

Total velocities relative to the centroid of the GNSS station triangle (mm/yr), with the translation vector removed

Site 1 azimuth and velocity relative to the triangle centroid (mm/yr)

{site1,(site1NormAzNSpeed[[1]]+centroidGridConvergence),site1NormAzNSpeed[[2]]}

Site 2 azimuth and velocity relative to the triangle centroid (mm/yr)

{site2,(site2NormAzNSpeed[[1]]+centroidGridConvergence),site2NormAzNSpeed[[2]]}

Site 3 azimuth and velocity relative to the triangle centroid (mm/yr)

{site3,(site3NormAzNSpeed[[1]]+centroidGridConvergence),site3NormAzNSpeed[[2]]}

Rotation

Rotation angular velocity and uncertainty (nano-radian/yr)

N[angVelNanoRadYr]

Rotation rate uncertainty (nano-radian/yr)

N[uncertRot]

Direction of rotation

If[angVelNanoRadYr<0,"clockwise","counter-clockwise"]

Lagrangian strain tensor

Instantaneous Lagrangian strain tensor ((Subscript[ε, ij]))

N[MatrixForm[m6]]

Uncertainties for instantaneous Lagrangian strain tensor ((Subscript[ε, ij]))

uncertainties={{uncertExx,uncertExy},{uncertExy,uncertEyy}};

N[MatrixForm[uncertainties]]

2D strain ellipse

Maximum extension (nano-strain)

N[e1/(10^(-9))]

Azimuth of maximum extension (degrees)

N[{azimuthS1HAxis1,azimuthS1HAxis2}]

Minimum extension (nano-strain)

N[e2/(10^(-9))]

Azimuth of minimum extension (degrees)

N[{azimuthS2HAxis1,azimuthS2HAxis2}]

Shear strain

Maximum shear strain (nano-strain)

N[maxShearStrain/(10^(-9))]

Area strain

Area strain (nano-strain)

N[areaStrain/(10^(-9))]

Invariants (nano-strain)

First invariant

N[firstInvariant/(10^(-9))]

Second invariant

N[secondInvariant/(10^(-9))]

Third invariant

N[thirdInvariant/(10^(-9))]

Plot of strain ellipse (exaggerated by 1 million times)

Plot of the infinitesimal horizontal strain ellipse, exaggerated by 1 million times.  The S1h axis is the red line, and S2h is blue.

summaryPlot

References

Allmendinger, R.W., Cardozo, N, and Fisher, D.M., 2012, Structural geology algorithms, vectors and tensors : Cambridge University Press, 289 p., ISBN 978 - 1 - 107 - 40138 - 9.

Cardozo, N., and Allmendinger, R.W., 2009, SSPX-- A program to compute strain from displacement/velocity data : Computers and Geosciences, v.35, p.1343 - 1357, doi : 10.1016/j.cageo .2008 .05 .008.

Means.W.D., 1976, Stress and strain - Basic concepts of continuum mechanics for geologists : Englewood Cliffs, New Jersey, Prentice - Hall, Inc., 339 p.

Savage, J.C., Gan, W., and Svarc, J.L., 2001, Strain accumulation and rotation in the Eastern California Shear Zone : Journal of Geophysical Research, v.106, B10, P.21, 995 - 22, 007.

Resources

Information about UNAVCO is available online via http : // www.unavco.org/

The full public data holdings of UNAVCO are available via their "Data Archive Interface Version 2" at http : // www.unavco.org/data/gps - gnss/data - access - methods/dai2/app/dai2.html.Information about the GPS - GNSS data holdings in particular is available via http : // www.unavco.org/data/gps - gnss/gps - gnss.html

General information about reference frames can be found at https : // www.unavco.org/software/visualization/GPS - Velocity - Viewer/GPS - Velocity - Viewer.html

The Universal Transverse Mercator coordinate system is described at https : // en.wikipedia.org/wiki/Universal_Transverse _Mercator _coordinate _system and https : // pubs.usgs.gov/fs/2001/0077/report.pdf and https : // gisgeography.com/utm - universal - transverse - mercator - projection/among other references.

The National Geodetic Survey' s NGS Coordinate Conversion and Transformation Tool (NCAT) is available via https : // www.ngs.noaa.gov/NCAT/

An online eigensystem calculator is available from Wolfram | Alpha at http : // www.wolframalpha.com/entities/calculators/eigenvalue_calculator/gv/9 h/bi/

A Wolfram | Alpha widget for evaluating a number of characteristics of a given matrix, including its eigensystem, is available online at http : // www.wolframalpha.com/widgets/view.jsp?id = ca851e9f71b25d2d588b7e0d4de3d627

A stand - alone Macintosh application called SSPX performs this same analysis in a more comprehensive manner for 3 or more GPS sites, including estimates of uncertainty.SSPX is available for academic or research use via the Mac App Store (https : // itunes.apple.com/us/app/sspx/id850389108); also visit Nestor Cardozo' s website at http : // www.ux.uis.no/~nestor/work/programs.html.