9. Coordinates, Projections, and Grids¶
9.1. Synopisis¶
- Review of coordinate systems and projections of the sphere onto the 2d placne
- Discussion about lengths and areas in finite volume grids used by ESMs
[1]:
%run -i figure-scripts/init.py
9.2. Coordinate systems¶
A coordinate system allows us to (uniquely) specify points in space. You should be familiar with the Cartesian 2d coordinate system where, by convention, (x,y) are the normal distances, measured in the same units, from two perpendicular lines through the origin.
We live in a 3d world (referring to spatial dimensions) and so require 3 numerical values to label a point in space. In Cartesian coordinates they might be (x,y,z) referenced to the center of the Earth, with z being height above the equatorial plane (positive in the direction of the North pole), y being the distance from a plane through the poles and a reference meridian, and x the distance from the plane perpendicular to both other planes. The equations of motion used by many models are often derived starting in this Cartesian coordinate system. However, these Cartesian coordinates are inconvenient to use in practice because we live on the surface of a sphere and “up”, as defined by gravity, is sometimes increasing z (at the North Pole), and sometimes changing x or y (at the Equator).
9.2.1. Spherical coordinates¶
In ESMs we typically use spherical coordinates, \(\lambda\), \(\phi\) and \(r\), where \(\lambda\) is “longitude”, a rotation angle eastward around the poles starting at a reference meridian; \(\phi\) is “latitude”, an elevation angle from the Equatorial plane (positive in Northern hemisphere), and \(r\) is the radial distance from the center of the Earth. \(\lambda,\phi,r\) are related to Cartesian \(x,y,z\) by some simple relations:
Note that \(r, x, y, z\) are all in the same units (eg. kilometets or meters) and \(\lambda,\phi\) are angles usually given in degrees or radians.
Coordinate singularities: At the North and South poles of the coordinate system, \(\phi = \pm \pi/2 = \pm 90^\circ\), all values of longitude refer to the same point. There is no “east” when you are positioned at the pole. This has many consequences, but one of the more fundamental is that spherical coordinates are not a good coordinate to use to design a discretization of the spherical domain.
Periodic coordinates: While a tuple of longitude, latitude and radius unambiguously define a point in space, given a point in space you there are multiple valid longitudes that refer to the same point. Longitude is cyclic (\(\pm360^\circ\) is equivalent to \(0^\circ\)). This can cause problems in practice, particularly for plotting spherical data for which effort is sometimes needed to handle the periodicity.
9.2.2. Geographic coordinates¶
We live on the surface of the Earth and to precisely refer to points near the Earth’s surface requires a properly defined geographic coordinate system. A common choice of coordinates is latitude, longitude and altitude, where altitude is height above a particular surface. Unfortunately the Earth is not spherical and that reference surface is better approximated as an ellipsoidal.
In order to be unambiguous about the definition of coordinates, map-makers choose a reference ellipse with a agreed upon scale and orientation. They then choose the most appropriate mapping of the spherical coordinate system onto that ellipsoid, called a geodetic datum. A widely used global datum includes the World Geodetic System (WGS 84), the default datum used for the Global Positioning System. When you are given a latitude-longitude pair of values, strictly speaking without the geodetic datum, the is some ambiguity about the actual physical point being referred to. For ESMs, the datum is rarely provided and this is because ESMs almost universally approximate the Earth as a sphere and use spherical coordinates for referencing locations. This means some approximation is required when comparing real-world positions and model positions.
The latitude and longitude using these horizontal datums are the spherical coordinates of the point on an ellipse. It you draw a straight line from the point on the ellipsoidal to the center, if passes through all spheres co-centered with the same latitude and longitude.
Different datum have different reference points and scales, and so longitude and latitude can differ between geodetic datum.
9.3. Projections¶
To view data covering the surface of a sphere, or the Earth, we have to project that 3d surface into 2d. Imagine peeling the rind off an orange in one piece and then trying to flatten it onto a table top; the curvature in the peel requires you to distort the rind or make cuts, in order to flatten it fully. This is the function of the map projections and distortion is inevitable. Some projection preserve properties such as relative angles between lines, or relative area, but there is no projection of the surface of the sphere that can avoid distortion of some form.
A projection maps the longitude and latitude of spherical coordinates into a new coordinate system. Very confusingly, sometimes the projection coordinates will be called longitude and latitude too! The projection coordinates are meaningless unless you know what the projection is so you often find a reference to the projection in the meta-data of coordinates; it means the longitude and latitude are not spherical coordinate but projection coordinates.
[2]:
%run -i figure-scripts/some-projections.py
Figure 1: The colored circles in these plots are circles in the tangent plane to the sphere and projected onto the surface. The various projections can distort the circles. The circles are separated zonally or meridionally by 60\(^\circ\). In 1a, a perspective image of the sphere, the circles appear non-circular because of the viewing angle. The blue circle appears circular because we are viewing it from directly overhead. The projection in 1b is the easy to use Plate-Carrée projection, a “lat-lon” plot, in which circles are stretched zonally with distance from the equator. 1d shows the Mercator projection in which circles remain circles but are expanded in size away from th equator. 1c shows the Robinson projection which compromises between the two. The purple dashed lines is a straight line in latitude-longitude coordinates, and the yellow dashed line is a straight line in the Mercator coordinates. The cyan dashed line is a great arc, and is straight in the perspective view because we are viewing it from directly overhead.
The two most useful projections are the equirectangular and Mercator projections.
9.3.1. Equirectangular projection¶
This is the simplest projection, sometimes thought of a non-projection which is incorrect. In general it takes the form
The origin of the plot, \((x,y)=(0,0)\) corresponds to \((\lambda,\phi)=(\lambda_0,\phi_0)\). The \(\cos \phi_0\) term is a constant and the most common choice of \(\phi_0=0\) gives the plate carrée projection, which means “flat square” in French. In this case, the projection is simply
Distances in the y-direction are proportional to the meridional direction on the sphere, but the x-direction are stretched zonally, more so further from the equator. This is apparent in the orange and green circles in figure 1b, where the heights or the loops are the same as circles on the equator but the width is markedly increased.
In the cartopy package, this projection is called “Plate-Carrée” which is French for flat square. Other names for this projection are equidistant cylindrical projection and geographic projection. See https://en.wikipedia.org/wiki/Equirectangular_projection.
9.3.2. Mercator projection¶
The Mercator projection has the same stretching in the x-direction as the equirectangular projection but, in order to preserve shape, it also stretches the y direction so that infinitesimal elements are stretched isotropically (the y-stretching is equal to the x-stretching).
At the polar singularities, the x-stretching is infinite so y becomes infinite and the Mercator projection can never show the poles. See https://en.wikipedia.org/wiki/Mercator_projection.
9.3.3. Lines¶
A length of a line between two points is a function of the path taken. On the surface of sphere, the shortest path between two given points is a great arc. A great arc does not appear straight in many projections. Unfortunately, many grid calculations use a great arc for the length of a line between nodes on a model grid, which can be inconsistent with the constraints or assumptions about the grid.
The dashed curves in figure 1 are “straight” lines between two points in various projections. The cyan dashed curve is a great arc. The purple dashed curve is a straight line in the Plate-Carree projection (latitude-longitude space) and the yellow dashed curve is a straight line in the Mercator projection. All are curved in most other projections. To describe a straight line in some projection then the projection must be known, irrespective of the coordinate system defining the end points. That is, we can define the end points of the line in latitude-longitude coordinates but say a line is straight in the Mercator projection, and by so doing unambiguously define that line.
In the Mercator projection, the length of a line is \(\frac{R}{\cos \alpha} \Delta \phi\) where \(\tan \alpha = \frac{\Delta y}{\Delta x}\)
9.4. ESM grids¶
Many ESMs use quadrilateral grids to discretize the surface of the sphere. The following discussion also applies to fully unstructured grids built from polygons but here we use quadrilateral grids for simplicity. There are also grids that have cuts and joins but here we’ll stick to space-filling grids that are logically rectangular, meaning they can be stored in rectangular arrays in computer memory and referenced with a pair of indices (\(i,j\) by convention).
A quadrilateral grid is a rectangular mesh of adjacent quadrilateral cells that share edges and vertexes. Although the mesh and the cell are logically rectangular they might be physically curvilinear. From the grid we require positions of nodes, distances along edges, and areas of cells.
If we choose a coordinate system with which to record the locations of mesh nodes, say spherical latitude-longitude with appropriate definitions, then we can unambiguously define those node locations. We could describe the exact same grid using a different coordinate system, say 3D Cartesian coordinates. The physical positions of the nodes of the grids are part of what define the grid, but the choice of coordinates with which we describe those positions does not change the grid.
The edges of each cell are a curve between two adjacent nodes but the particular path of the curve has to be defined. Different paths will have different lengths. Similarly, the particular paths of the cell edges will determine the cell area. Thus the path of the cell edges is a fundamental component of a model grid needed for calculating the lengths and areas on a grid.
9.4.1. Simple spherical coordinate grid¶
Before we discuss the best choice for defining a curve between points, let’s briefly define a simple spherical-coordinate grid. The mesh is formed of lines of constant longitude and lines of constant latitude.
Let \(i \in 0,1,\ldots, n_i\) and \(j \in 0,1,\ldots, n_j\), then node \(i,j\) is at longitude \(\lambda_i\) and latitude \(\phi_j\) where \(\lambda_i=\lambda_0 + i \Delta \lambda\), \(\phi_j=\phi_0 + j \Delta \phi\).
Here, \(\Delta \lambda\) and \(\Delta \phi\) are grid spacings. In practice, these can be smooth functions of \(i\) and \(j\) respectively but here we treat them as constant.
An example simple spherical grid is shown below. The red dots are the nodes of the mesh with positions \(\lambda_i,\phi_j\). The dashed lines are the cell edges that for a regular net. Notice that in the Plate-Carrée projection the grid is regular because the grid-spacing in constant in longitude-latitude coordinates.
The lengths and areas of the grid are measured on the surface of sphere. We defined the edges to be either lines of constant longitude or latitude. Using spherical geometry, the length of a meridionally-oriented (constant longitude) cell edge is \(R \Delta \phi\). For a zonally-oriented edge at constant latitude \(\phi_j\), the length is \(R \Delta \lambda \cos \phi_j\). The area of a cell labelled \(i+\frac{1}{2},j+\frac{1}{2}\) bounded by four edges is \(R^2 \Delta \lambda \left( \sin \phi_{j+1} - \sin \phi_j \right)\).
The metric factors for this grid are the same as for a Plate-Carrée projection because we are defining the paths of the cell edges to be straight in the Plate-Carrée projection. The use of the Plate-Carrée coordinates for position, namely longitude and latitude, is a happy coincidence which means everything, positions and metrics, are defined by one projection.
[3]:
%run -i figure-scripts/simple-spherical-grid.py
Figure 2: A simple spherical grid view in different projections. The red dots are the nodes of the mesh. The dashed lines are the grid of cell edges. Here the grid spacing is constant in latitude-longitude coordinates which is why the grid looks regular in the Plate-Carrée projection.
In the above example of a simple spherical grid we used spherical geometry combined with the specification of the paths of cell edges to calculate everything. If we were given just the positions of nodes a spherical grid, we could re-calculate the grid lengths and areas given the knowledge that the cell edges are straight in the Plate-Carrée projection.
9.4.2. Not to use great arcs¶
If is quite common for the projection for the cell edges to be omitted from a grid file. Sometimes, you will encounter some software that was written to handle arbitrary grids that does not know about the projections. In this instance, developers often make a choice for how to calculate a length between two points. There are two immediately simple choices we can make: i) the shortest curve, or ii) choose a projection in which to draw a straight line. There are other ways to choose and define the curve forming edges but they are not as useful for out purposes.
As discussed above, the shortest line is a great arc but it turns out that choosing great arcs does not result in orthogonal meshes without a particular distributions of nodes. For example, if you consider the example spherical grid defined above, using a great arc to approximate the distance along a latitude circle will underestimate the distance. If you were to then calculate the actual paths implied by the great arcs, the grid is not orthogonal at the nodes.
If no projection is given, then straight lines in the Plate-Carrée projection make a reasonable approximation that at least works for Plate-Carrée and Mercator. Many grids a composed of patches of grids that used different projections but have to match at the joins and it seems the Plate-Carrée is often the common denominator.
9.5. Contents of grid files¶
As mentioned earlier, because most ESM treat the Earth as spherical, they rarely bother to choose a datum for the geographic coordinate system. Further, the longitudes and latitudes may not be spherical coordinates but instead may be the result of a projection.
The habit of not specifying the projection for cell edges is somewhat alleviated by the general pattern of providing the numerical values for lengths and areas as part of the grid specification. For example, see https://www.researchgate.net/publication/228641121_Gridspec_A_standard_for_the_description_of_grids_used_in_Earth_System_models. For many calculations, including integrals and gradients of scalars, having the numerical values of lengths and areas is sufficient. For vector operations, the angle of grid-lines is needed.
The CF convention requires the grid to be provided in all files. That convention supports specification of mappings (projections) but does not distinguish the projection for paths from the projections for coordinates. by default, longitude and latitude are the true geographic coordinates. And the required grid information is limited to the node positions. Further the nodes correspond to the data locations which is a finite-difference perspective of a mesh, and quite different from the finite volume perspective that most ESMs assume.