2a:Tf58a,heightstxt = md`## Understanding the geoid and "height" The [*geoid*](https://oceanservice.noaa.gov/facts/geoid.html) is one approximation of Earth's shape due to gravity, roughly coinciding with mean sea level. Because the inside of our planet is non-uniform, gravity is not constant at the surface. If Earth were completely covered in an ocean, the water would swell over gravity highs (where there is more mass below the surface) and depress in areas of gravity lows (where there is less mass below the surface). The geoid approximates this hypothetical ocean surface, representing a geopotential surface. To think of it another way, if we were able to walk on the geoid, we would experience a constant value of gravity. The height of the geoid is often reported relative to the internationally recognized [WGS84 *reference ellipsoid*](https://geodesy.noaa.gov/datums/index.shtml). The GPS in our phones, for example, records elevations relative to this ellispoid. Sometimes, it is useful to know elevations relative to mean sea level—this is where the geoid height calculation comes in: to convert ellipsoidal heights to heights above the geoid (aka *orthometric heights*). The NGA ([National Geospatial-Intelligence Agency](https://earth-info.nga.mil/index.php?action=home)) publishes models that describe the earth's gravitational field (see the Earth Gravitational Models [here](https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84)), and updates them as we learn more through satellite geodesy. It is becoming standard to use EGM2008 over its predecessors, but a 2020 model is soon to become standard. Here is a helpful sketch of what these surfaces might look like relative to Earth's topography:

drawing

*Image from the [National Center for Polar and Ice Research](https://ncpor.res.in/pages/researchview/11)*

**Terminology definitions:** - **Reference Ellipsoid:** Here, the reference ellipsoid refers to the datum surface derived from the WGS84 reference system. - **Geoid:** The gravitational equipotential surface that approximates global mean sea level, given as the height above the reference ellipsoid. - **Geoid Undulation:** Height of the geoid above the reference ellipsoid. - **Geoid Height**: Same as geoid undulation. - **Orthometric Height:** Height above the geoid. `2b:T75e,calculationtxt = md`## Understanding the geoid calculation Computing the height of a geoid from a set of potential coefficients provided in an Earth Gravitational Model (EGM), as done here, is a non-unique problem, dependent on several assumptions. If you'd like to know more about the calculation, a core paper that outlines the methodology can be found here: [*Use of Potential Coefficient Models for Geoid Undulation Determination Using a Spherical Harmonic Representation of the Height Anomaly/Geoid Undulation Difference*](https://link.springer.com/article/10.1007%252Fs001900050096,) by R.H. Rapp, Journal of Geodesy, 1996. Further discussion of the assumptions in computing a geoid from an EGM is presented in [*There is no such thing as "The" EGM96 geoid: Subtle points on the use of a global geopotential model*](https://www.ngs.noaa.gov/PUBS_LIB/EGM96_GEOID_PAPER/egm96_geoid_paper.html), by Dru A. Smith, 1998. A good description of errors in the geoid height calculations, including those introduced by the cubic interpolation method, can be found on the [GeographicLib](https://geographiclib.sourceforge.io/html/geoid.html#geoidinterp) page. In sum, the errors from the calculations presented here (when compared against the NGA geoid [test datasets](https://observablehq.com/@bmershon/geoid-height#tests)) stem from (1) the method of interpolation, (2) the height ambiguity resulting from the fact that the geoid and ellipsoid are not parallel surfaces, and (3) the grid format of the geoid harmonic coefficients, which require quantization of height to 3 mm. A comparison of errors from this calculator with the NGA test datasets is performed in a companion [Notebook](https://observablehq.com/@bmershon/geoid-height) created by Brooks Mershon and Cassie Hanagan at UNAVCO. An associated map highlights areas where the geoid height calculations differ the most. `2c:T574,origintxt = md`## Code origin The original UNAVCO/NSF GAGE tool ran on a NGA ([National Geospatial-Intelligence Agency](https://earth-info.nga.mil/index.php?action=home)) fortran code for the geoid height calculation. The geoid height calculation was recoded in C++ as part of a larger library of geodetic functions in [GeographicLib](https://geographiclib.sourceforge.io/). The [GeographicLib::Geoid](https://geographiclib.sourceforge.io/html/classGeographicLib_1_1Geoid.html) class and associated [Geoid.cpp](https://geographiclib.sourceforge.io/html/Geoid_8cpp.html) were translated to JavaScript by [Kim Vandry](vandry@TZoNE.ORG), available [here](https://github.com/vandry/geoidheight). All products are under the [MIT License](https://opensource.org/licenses/MIT). This notebook is run with core functions from another [Notebook](https://observablehq.com/@bmershon/geoid-height) created by Brooks Mershon and Cassie Hanagan at the GAGE Facility, which makes Kim Vandry's JavaScript port of GeographicLib::Geoid available in Observable Notebooks like this one. That notebook includes a [suite of test cases](https://observablehq.com/@bmershon/geoid-height#tests) and further documentation to validate the accuracy of this calculator. `2d:[{"id":"6e95cb17aaefc544c9c7d7fe8afc92cd79df09fdf32bd2788cb7075f2baa1da257ccf9091b239176fbe2b62d926f76227883f55b0649a1657969b1fc9b30a20f","href":"https://static.observableusercontent.com/files/6e95cb17aaefc544c9c7d7fe8afc92cd79df09fdf32bd2788cb7075f2baa1da257ccf9091b239176fbe2b62d926f76227883f55b0649a1657969b1fc9b30a20f","name":"unav_50_real-time_stations@1.csv","type":"text/csv","lastModified":1624985858699.833,"size":237},{"id":"0ca9148938fe656b24d3cd9474fc7190de67f0c2fa0c6940f93f874da1572341e3691665fd40b708cdd5573317e4b18f96bfc1d3617464cb0dc2d8b607c3bddd","href":"https://static.observableusercontent.com/files/0ca9148938fe656b24d3cd9474fc7190de67f0c2fa0c6940f93f874da1572341e3691665fd40b708cdd5573317e4b18f96bfc1d3617464cb0dc2d8b607c3bddd","name":"world-atlas-110m.json","type":"application/json","lastModified":1624299310527.743,"size":102116}] 29:{"title":"Geoid Height Calculator","cells":[{"id":0,"value":"introtxt = md`# Geoid Height Calculator\n\n© [NSF GAGE Facility](https://www.unavco.org/) by Cassie Hanagan & Brooks Mershon\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2578,"value":"**See also**: [Earth Gravitation Model](https://observablehq.com/@earthscope/egm) for more visual explanation.","pinned":false,"hidden":false,"mode":"md"},{"id":4,"value":"overviewtxt = md`## Overview\n\nThis tool computes the geoid height and orthometric height of a point (or points). \n\nUser-specified inputs include:\n 1. The [Earth Gravitational Model](https://earth-info.nga.mil/index.php?dir=wgs84&action=wgs84) to compute the geoid. The number following EGM represents the model year, the number following the dash denotes the model grid spacing in minutes of a degree.\n 2. The option for cubic interpolation between grid points. Otherwise no interpolation is performed. *See the [GeographicLib](https://geographiclib.sourceforge.io/html/geoid.html#geoidinterp) information pages for more info on the interpolation method.* \n 3. The longitude, latitude, and height of a point(s) above the WGS84 reference ellipsoid.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1262,"value":"calc = md`## Calculator: Input and output\n\n**Select the geoid model to use and choose whether you want cubic interpolation between grid points:**`","pinned":false,"hidden":false,"mode":"ojs"},{"id":345,"value":"viewof modelSelection = Inputs.select(models, {\n value: models[0],\n label: html`1. Geoid Model`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":331,"value":"viewof cubic = Inputs.toggle({\n label: html`2. Cubic interpolation`,\n value: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1145,"value":"instrtxt = md`**Import CSV or copy & paste point data**\n\n(1) Upload a comma-separated file with longitude (decimal degrees), latitude (decimal degrees), and ellipsoid height (m) column headers. \n\n**or**\n\n(2) Copy and paste comma-separated longitude (decimal degrees), latitude (decimal degrees), and ellipsoid height (m) values without headers.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1065,"value":"viewof csvData = csvInput({\n initialUrl: await FileAttachment(\"unav_50_real-time_stations@1.csv\").url(),\n label: \"CSV File:\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2764,"value":"note3 = md`Default coordinates below include:\n- Boulder, CO\n- Mount Everest\n- a gravity hole in the Indian Ocean.\n\n*Edit and submit new values below as desired.*`","pinned":false,"hidden":false,"mode":"ojs"},{"id":513,"value":"viewof coordinatesInput = textarea({\n title: \"Coordinates (Longitude, Latitude, Ellipsoid Height)\",\n value: defaultCsv, // Random sample of points and elevation\n rows: 10,\n cols: 60,\n submit: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":480,"value":"earth = Plot.plot({\n title: \"Your Coordinates\",\n color: {\n legend: true,\n type: \"diverging\",\n domain: [100, -100]\n },\n width: width,\n projection: { type: \"equal-earth\" },\n style: \"overflow: visible;\", // allow dots to escape\n marks: [\n Plot.geo(land, {\n fill: \"currentColor\",\n stroke: \"currentCOlor\",\n fillOpacity: 0.2\n }),\n Plot.raster(matrix, {\n clip: \"sphere\",\n fill: \"MAG_IGRF90\",\n interpolate: \"barycentric\",\n fill: (d, i) => d[2],\n fillOpacity: opacityEarth\n }),\n Plot.graticule(),\n Plot.dot(coords, {\n clip: \"sphere\",\n opacity: opacityEarth,\n r: 4,\n stroke: \"yellow\",\n fill: \"black\"\n }),\n Plot.tip(\n matrix,\n Plot.pointer({\n x: (d) => d[0],\n y: (d) => d[1],\n title: (d) => `lon/lat (${d[0]}, ${d[1]})`\n })\n )\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2651,"value":"viewof opacityEarth = Inputs.range([0, 1], {\n label: \"Opacity\",\n value: 0.7,\n step: 0.1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2803,"value":"tableExplanation = md`The table below contains your given ellipsoid height samples and their corresponding geoid and orthometric heights.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1678,"value":"viewof solutionTable = Inputs.table(csvOut, {\n format: {\n Longitude: (x) => html`${d3.format(\"D\")(x)}`,\n Latitude: (x) => html`${d3.format(\"D\")(x)}`,\n \"Ellps Height (m)\": (x) => html`${d3.format(\"D\")(x)}`,\n \"Geoid Height (m)\": (x) =>\n html`${x}`,\n \"Ortho Height (m)\": (x) =>\n html`${x}`\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1244,"value":"savemetxt = md`Click the button below to save and download a zip file with the table in CSV format and a README describing the calculation input parameters.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1922,"value":"downloadzip = download(\n zip,\n 'GeoidHeight_' + d3.timeFormat(\"%Y-%m-%d\")(new Date())\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"$2a","pinned":false,"hidden":false,"mode":"ojs"},{"id":2581,"value":"## Visualizing Geoid Undulation","pinned":false,"hidden":false,"mode":"md"},{"id":2612,"value":"_[The Earth Gravitation Model Notebook](https://observablehq.com/@earthscope/egm) can help us visualize the geoid undulation compared to familiar land masses and oceans. The globe below defaults to a gravity hole in the Indian Ocean as a point of interest._","pinned":false,"hidden":false,"mode":"md"},{"id":2586,"value":"import { globe, longitude, latitude, matrix, meridianPlot, parallelPlot } with {\n selectedLongitude as longitude,\n selectedLatitude as latitude,\n selectedPojection as projection,\n plotHeight as plotHeight\n} from \"@earthscope/egm\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2595,"value":"viewof selectedLongitude = Inputs.range([-180, 180], {\n label: \"longitude\",\n step: 1,\n value: 79\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2602,"value":"viewof selectedLatitude = Inputs.range([-90, 90], {\n label: \"latitude\",\n step: 1,\n value: 5\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2607,"value":"viewof selectedPojection = Inputs.select([\"Orthographic\", \"Equal-Earth\"], {\n label: \"Projection\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2769,"value":"The rasterized (interpolated) values for each latitude and longitude sample on the globe below show the relatively smooth undulation of the geoid as a visualization of extended mean sea level across the entire globe.","pinned":false,"hidden":false,"mode":"md"},{"id":2589,"value":"globe","pinned":false,"hidden":false,"mode":"ojs"},{"id":2835,"value":"meridianPlot","pinned":false,"hidden":false,"mode":"ojs"},{"id":2841,"value":"parallelPlot","pinned":false,"hidden":false,"mode":"ojs"},{"id":407,"value":"$2b","pinned":false,"hidden":false,"mode":"ojs"},{"id":297,"value":"$2c","pinned":false,"hidden":false,"mode":"ojs"},{"id":372,"value":"reftxt = md`## References \n\nC. F. F. Karney, GeographicLib, Version 1.51 (2019-09-24), https://geographiclib.sourceforge.io/1.51\n\nRapp, R. H. (1997). Use of potential coefficient models for geoid undulation determinations using a spherical harmonic representation of the height anomaly/geoid undulation difference. *Journal of Geodesy, 71(5)*, 282-289.\n\nSmith, D. A. (1998). There is no such thing as \"The\" EGM96 geoid: Subtle points on the use of a global geopotential model. *IGeS Bulletin N. 8*, Int. Geoid Service, 17-27.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":310,"value":"md`### Appendix`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2868,"value":"pageview","pinned":false,"hidden":false,"mode":"ojs"},{"id":2847,"value":"plotHeight = 200","pinned":false,"hidden":false,"mode":"ojs"},{"id":357,"value":"md`** I. Code Imports**`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2496,"value":"import { Copier } from \"@mbostock/copier\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":311,"value":"import { models, geoid } from \"@unavco/geoid-height-function/2\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":518,"value":"import { select, textarea, slider, color } from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1059,"value":"import { csvInput } from \"@unavco/file-input-with-csv-support\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1241,"value":"import { download } from \"@mbostock/lazy-download\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2344,"value":"import { legend, swatches } from \"@d3/color-legend\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1271,"value":"md` ** II. Non-geoid calculation variables and functions**\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2559,"value":"function toc(selector = \"h2\", heading = \"Jump to:\") {\n return Generators.observe((notify) => {\n let headings = [];\n\n function observed() {\n const h = Array.from(document.querySelectorAll(selector));\n if (h.length !== headings.length || h.some((h, i) => headings[i] !== h)) {\n notify(\n html`${heading} ${Array.from((headings = h), (h) => {\n return Object.assign(\n html`${DOM.text(\n h.textContent\n )} `,\n { onclick: (e) => (e.preventDefault(), h.scrollIntoView()) }\n );\n })}`\n );\n }\n }\n\n const observer = new MutationObserver(observed);\n observer.observe(document.body, { childList: true, subtree: true });\n observed();\n return () => observer.disconnect();\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":491,"value":"graticule = d3.geoGraticule10()","pinned":false,"hidden":false,"mode":"ojs"},{"id":575,"value":"project = {\n return d3[\"geoNaturalEarth1\"]()\n .rotate([0, 0])\n .precision(0.1)\n .fitSize([width, width / 1.8], { type: \"Sphere\" });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":494,"value":"land = topojson.feature(world, world.objects.land)","pinned":false,"hidden":false,"mode":"ojs"},{"id":496,"value":"world = FileAttachment(\"world-atlas-110m.json\").json()","pinned":false,"hidden":false,"mode":"ojs"},{"id":498,"value":"topojson = require(\"topojson@3\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":499,"value":"d3 = require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":593,"value":"countries = topojson.feature(world, world.objects.countries)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1697,"value":"function rowsToObjects(rows, headers) {\n return rows.map(r =>\n r.reduce((o, c, i) => Object.assign(o, { [headers[i]]: c }), {})\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1866,"value":"JSZip = require(\"jszip@3/dist/jszip.min.js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":2365,"value":"colorscale = d3.scaleDiverging(\n [d3.min(geoid_heights), 0, d3.max(geoid_heights)],\n d3.interpolateRdBu\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":380,"value":"md`**III. Geoid calculator variables and functions**\n\nThe user-specified inputs are run through the geoid function in the following cells.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1123,"value":"// lon, lat\ncoords = coordinatesInput\n .trim()\n .split(/\\n/)\n .map(\n line =>\n line\n .split(/,\\s*/)\n .slice(0, 2) // elevation discarded\n .map(n => +n) // number\n )","pinned":false,"hidden":false,"mode":"ojs"},{"id":1294,"value":"heights = coordinatesInput\n .trim()\n .split(/\\n/)\n .map(\n line =>\n line\n .split(/,\\s*/)\n .slice(-1) // lon, lat discarded\n .map(n => +n) // number\n )","pinned":false,"hidden":false,"mode":"ojs"},{"id":2635,"value":"defaultCsv = `-105.2705,40.0150,1621\n86.92528,27.98833,8848.86\n79,5,0`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2675,"value":"fetchElevation = function (lons, lats) {\n return fetch(\n `https://api.open-meteo.com/v1/elevation?latitude=${lats.join(\n \",\"\n )}&longitude=${lons.join(\",\")}`\n ).then((response) => response.json());\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1280,"value":"geoid_heights = {\n let gheight = [];\n for (var i = 0; i < coords.length; i++) {\n gheight.push(await new geoidHeight(coords[i][0], coords[i][1], cubic));\n }\n return gheight;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1278,"value":"ortho_heights = {\n let oheight = []; //; var oheight = [];\n for (var i = 0; i < coords.length; i++) {\n oheight.push(heights[i] - geoid_heights[i]);\n }\n return oheight;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2312,"value":"// reverse coordinates for mapping\nmap_coords = {\n let n = [];\n for (var i = 0; i < coords.length; i++) {\n n.push([coords[i][0], coords[i][1], geoid_heights[i]]);\n }\n return n.map(line => line.slice(0, 2).reverse());\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1527,"value":"// filter coordinates based on whether the geoid is above the ellipsoid\nabove_coords = {\n let n = [];\n for (var i = 0; i < coords.length; i++) {\n n.push([coords[i][0], coords[i][1], geoid_heights[i]]);\n }\n return n\n .filter(function(arr) {\n return arr[2] > 0;\n })\n .map(line => line.slice(0, 2).reverse());\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1435,"value":"below_coords = {\n let n = [];\n for (var i = 0; i < coords.length; i++) {\n n.push([coords[i][0], coords[i][1], geoid_heights[i]]);\n }\n return n\n .filter(function(arr) {\n return arr[2] <= 0;\n })\n .map(line => line.slice(0, 2).reverse());\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1650,"value":"// Maps the height, ortho height, and geoid height to the coordinates,\n// and turns the array into an object with column headers\ncsvOut = rowsToObjects(\n coordinatesInput\n .trim()\n .split(/\\n/)\n .map(\n (line) =>\n line\n .split(/,\\s*/)\n .slice(0, 2) // elevation discarded\n .map((n) => +n)\n .reverse() // number\n ) // Left with lat, lon. Next: add heights\n .map((line) => line.reverse())\n .map((d, i) => {\n (d[2] = heights[i][0]),\n (d[3] = Number(geoid_heights[i].toFixed(2))),\n (d[4] = Number(ortho_heights[i].toFixed(2)));\n return d;\n }),\n [\n \"Longitude\",\n \"Latitude\",\n \"Ellps Height (m)\",\n \"Geoid Height (m)\",\n \"Ortho Height (m)\"\n ]\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1869,"value":"function readmeText() {\n return `© ${d3.timeFormat(\"%Y\")(new Date())} UNAVCO\n\n Created with the Geoid Height Calculator at:\n\n https://observablehq.com/@cehanagan/geoid-height-calculator\n\n The file \"geoid_height_calculator_${d3.timeFormat(\"%Y-%m-%d\")(\n new Date()\n )}.csv\"\n contains the user-input latitudes, longitudes, and heights (above the WGS84 ellipsoid), \n and calculated geoid and orthometric heights above the ${modelSelection} geoid.\n\n Cubic interpolation was set to ${cubic}.\n\n Contact brooks.mershon@unavco.org or @bmershon on GitHub with questions, concerns, or 🐛 reports.\n \n `;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1880,"value":"zip = {\n const zip = new JSZip();\n const folder = zip.folder(\n 'GeoidHeight' + d3.timeFormat(\"%Y-%m-%d\")(new Date())\n );\n folder.file('README.txt', readmeText());\n folder.file(\n 'geoid_height_' + d3.timeFormat(\"%Y-%m-%d\")(new Date()) + '.csv',\n d3.csvFormat(solutionTable)\n );\n return zip.generateAsync({ type: 'blob' });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1015,"value":"function geoidHeight(lon, lat, cubic = false) {\n return geoid(modelSelection).then(\n // This function is called when the dataset is loaded (a promise resolves).\n function (calculate) {\n // The parameter to this function is a function that takes a lat/long and provides a\n // height using the loaded dataset.\n return calculate(\n lat, // latitude\n lon, // longitude (be careful about the order of the arguments!)\n cubic // true or false\n );\n }\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2672,"value":"import { pickRandom } from \"@earthscope/random\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":2865,"value":"import { pageview } from \"@earthscope/pageview\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W2d","stdlib":"1","titling":"auto"}