2f:[] 2d:{"title":"GNSS & position & regression, oh my!","cells":[{"id":0,"value":"# GNSS & position & regression, oh my!","pinned":false,"hidden":false,"mode":"md"},{"id":318,"value":"This notebook attempts to *interactively* illustrate the findings published by Blewitt and Lavalleé in 2001 regarding bias in estimation of site velocities.","pinned":false,"hidden":false,"mode":"md"},{"id":308,"value":"md`Blewitt, G., and Lavallée, D., Effect of annual signals on geodetic velocity, J. Geophys. Res., 107( B7), doi:10.1029/2001JB000570, 2002.\n
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":210,"value":"## Least Squares\nIf we naively fit a line using least squares to compute a velocity from daily positions for a given GNSS ground station, the resulting velocity estimate may be biased by un-modeled annual harmonics which are often present in GNSS timeseries.\n\nBelow we have a simulated a daily GNSS time-series that includes an annual harmonic with an amplitude of **${amplitude} cm**. The simulated time-series does not include a velocity. However, the resulting line-fit indicates a negative trend or velocity. ","pinned":false,"hidden":false,"mode":"md"},{"id":34,"value":"viewof amplitude = Inputs.range([0, 2], {\n label: \"Amplitude (cm)\",\n step: 0.1,\n value: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":520,"value":"We can also model a velocity in our time series. See if you can create both a positive and negative slope using the slider below.","pinned":false,"hidden":false,"mode":"md"},{"id":151,"value":"viewof velocity = Inputs.range([-1, 1], {\n label: \"Velocity [cm/year]\",\n step: 0.1,\n value: 0.2 // cm\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":262,"value":"chartNoNoise = Plot.plot({\n width: width,\n x: {\n label: \"Days\",\n grid: true\n },\n y: {\n domain: [-3, 3],\n label: \"Displacement Up (cm)\",\n grid: true\n },\n marks: [\n Plot.dot(position, { x: (d) => d[0], y: (d) => d[1] }),\n Plot.linearRegressionY(position, {\n x: (d) => d[0],\n y: (d) => d[1],\n stroke: \"orange\",\n ci: 0.95\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":271,"value":"## Noise\nOf course, our simulation should probably have a little noise. Here's a Normal distribution that assigns a random value for each position observation above. We can control the \"noisiness\" by changing the standard deviation (`sigma`).","pinned":false,"hidden":false,"mode":"md"},{"id":477,"value":"viewof sigma = Inputs.range([0, 2], {\n label: \"Sigma\",\n step: 0.01,\n value: 0.7\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":416,"value":"normal_noise_plot = Plot.plot({\n width: width,\n y: {\n domain: [-3, 3],\n grid: true\n },\n marks: [\n Plot.dot(randomNormal, { x: (d) => d[0], y: (d) => d[1] }),\n Plot.linearRegressionY(randomNormal, {\n x: (d) => d[0],\n y: (d) => d[1],\n stroke: \"blue\",\n ci: 0.95\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":275,"value":"We can add the noise above to our model data to obtain this series of position observations with a real drift of **${d3.format(\".2r\")(velocity * 10)} mm** per year.","pinned":false,"hidden":false,"mode":"md"},{"id":5,"value":"exported_plot = Plot.plot({\n width: width,\n x: {\n label: \"Time\",\n grid: true\n },\n y: {\n domain: [-3, 3],\n label: \"Displacement (cm)\",\n grid: true\n },\n marks: [\n Plot.dot(\n position.map((d, i) => [d[0], d[1] + randomNormal[i][1]]),\n { x: (d) => d[0], y: (d) => d[1] }\n ),\n Plot.linearRegressionY(position, {\n x: (d) => d[0],\n y: (d) => d[1],\n stroke: \"orange\",\n ci: 0.95\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":524,"value":"## Phase Shift\nTry changing the phase of the harmonic using the toggle below. Notice anything about the orange line's slope?","pinned":false,"hidden":false,"mode":"md"},{"id":66,"value":"viewof phase = Inputs.range([0, 365], {\n label: \"Phase Shift [days]\",\n step: 1,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":522,"value":"How about if you change the length of the times series, does the slope change?","pinned":false,"hidden":false,"mode":"md"},{"id":142,"value":"viewof years = Inputs.range([1, 10], {\n label: \"Duration [years]\",\n step: 0.5,\n value: 3\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":582,"value":"## Bias\nOur simulated observations above modeled by a real drift, an annual signal, and some noise gives us the following slope/intercept equation.","pinned":false,"hidden":false,"mode":"md"},{"id":561,"value":"tex.block`\ndisplacement = ${d3.format(\".3r\")(solution.a)} t ${\n solution.b > 0 ? \"+\" : \"-\"\n} ${d3.format(\".3r\")(Math.abs(solution.b))}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":592,"value":"... where ${tex`t`} is our day since observations began and displacement is in centimeters.","pinned":false,"hidden":false,"mode":"md"},{"id":597,"value":"tex.block`\n${d3.format(\".3r\")(solution.a)} (365) = ${d3.format(\".3r\")(solution.a * 365)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":630,"value":"Compare this derived slope from our noisy phase-shifted data spanning **${years} years** above that has a real velocity of **${velocity} cm/year**.","pinned":false,"hidden":false,"mode":"md"},{"id":645,"value":"error_plot = Plot.plot({\n width: width,\n x: {\n label: \"Time\",\n grid: true\n },\n y: {\n domain: [-3, 3],\n\n label: \"Displacement (cm)\",\n grid: true\n },\n marks: [\n Plot.dot(\n position.map((d, i) => [d[0], d[1] + randomNormal[i][1]]),\n { x: (d) => d[0], y: (d) => d[1], fill: \"#ccc\" }\n ),\n Plot.line([\n [0, 0],\n [years * 365, years * velocity]\n ]),\n Plot.linearRegressionY(position, {\n x: (d) => d[0],\n y: (d) => d[1],\n stroke: \"orange\",\n ci: 0.95\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":636,"value":"## Data Span\nTo reiterate the conclusion put forth by Blewitt and Lavallée:\n\n### *\"We recommend that 2.5 years be adopted as a standard minimum data span for velocity solutions intended for tectonic interpretation or reference frame production and that we be skeptical of geophysical interpretations of velocities derived using shorter data spans.\"*","pinned":false,"hidden":false,"mode":"md"},{"id":694,"value":"Can you think of a way to measure velocity that is less sensitive to the annual signal?","pinned":false,"hidden":false,"mode":"md"},{"id":696,"value":"## Further Reading\n\n- Blewitt, G., Kreemer, C., Hammond, W. C., and Gazeaux, J. (2016), MIDAS robust trend estimator for accurate GPS station velocities without step detection, J. Geophys. Res. Solid Earth, 121, 2054– 2068, doi:10.1002/2015JB012552.","pinned":false,"hidden":false,"mode":"md"},{"id":296,"value":"---\n# Appendix","pinned":false,"hidden":false,"mode":"md"},{"id":540,"value":"solution = statistics\n .regressionLinear()\n .x((d) => d[0])\n .y((d) => d[1])\n .domain(d3.extent(domain))(\n position.map((d, i) => [d[0], d[1] + randomNormal[i][1]])\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":167,"value":"domain = d3.range(0, years * 365, 1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":403,"value":"randomNormal = domain.map((d) => [d, d3.randomNormal(0, sigma)()])","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"// \"Daily\" position solutions from a GNSS station\nposition = domain.map((day, i) => [\n day,\n simulatedTSFunc(day),\n]);\n\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":509,"value":"createSimulatedTSFunc = function(velocity, phase, amplitude) {\n return function(day) { return (velocity/365) * day + amplitude * Math.sin(((2*Math.PI)/365)*(day + phase))}\n};","pinned":true,"hidden":false,"mode":"ojs"},{"id":508,"value":"simulatedTSFunc = createSimulatedTSFunc(velocity, phase, amplitude);","pinned":true,"hidden":false,"mode":"ojs"},{"id":542,"value":"statistics = require(\"d3-regression\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W2f","stdlib":"1","titling":"auto"}