This post came about from thinking about how to properly project WRF output for visualisation using Cartopy. It is also relevant for anyone wondering how to properly interpolate to real world locations if high accuracy is needed. A follow-up post will deal with using WRF and Cartopy
Location, location, location
As with any earth-science visualisation, the most confusing part is understanding coordinate systems and projections. It should be simple – we have rectangular arrays representing variables and we wish to display them on a rectangular part of the screen. How does it get so complicated? When you set up a WRF simulation, you define a projection. That is, you specify how a part of the curved earth is going to be represented in rectangular arrays in memory and disk. A good article on projections can be found here. An important point to note is that WRF and all other climate models assume a spherical earth, which makes their code much simpler.
Latitude is ambigouous
The latitude of a point on the surface of a sphere is unambiguously the angle from the equator to that point. When you start considering the Earth as a ellipse or something lumpier …
- sphere, one parameter e.g. WRF R=6 370 000m
- ellipsoid, two paramaters e.g. WGS84 a=6 378 137m, b=6 356 752.3142m
- geoid, lumpy
… then it gives rise to some abiguities when definining latitude, shown schematically below image from here. According to Wikipedia “The terminology for latitude must be made more precise by distinguishing:
Spherical latitude: the angle between the normal to a spherical reference surface and the equatorial plane.
Geodetic latitude: the angle between the normal and the equatorial plane. The standard notation in English publications is φ. This is the definition assumed when the word latitude is used without qualification. The definition must be accompanied with a specification of the ellipsoid.
Geocentric latitude: the angle between the radius (from centre to the point on the surface) and the equatorial plane. There is no standard notation: examples from various texts include ψ, q, φ’, φc, φg.
Atsronomic latitude the angle between the equatorial plane and the true vertical (plumb line) at a point on the surface.
Geographic latitude must be used with care. Some authors use it as a synonym for geodetic latitude whilst others use it as an alternative to the astronomical latitude.”
WRF Preprocessing Georeferenced Data
An important but not well documented point is how georeferenced input data is handled in WRF. This includes terrain height, land-use data, and boundary conditions from GFS or other models. According to this post on using your own high-resolution data, “WPS considers the Earth to be a perfect sphere … default land-use data used by WPS are projected using the WGS_1984 coordinate system. Due to this reason, it is preferred to convert all the land data to the following projection GCS_WGS_1984”
From this, I infer:
- WPS treats all latitudes and longitudes as spherical, indentifying points on a sphere.
- Static data like land-use and terrain data are most likely georefrenced against the WGS84 ellipsoid.
- An implicit transformation is occuring in WPS, whereby the WGS elliposid is projected onto a sperical surface.
Something like this thanks, with each P –> P’
This is fine, as long as it is consistent between all input data. Once concern is the treatment of GFS or other global model data. How are the observations handled in the GFS data assimilation system? In the best case, the same approach is used, and we can sleep easily at night. What if a different transformation is being applied to transform observations to the spherical system used in GFS?
In that case, when we overlay meteorological data from GFS over the static data from WPS, there could be a displacement. Does it matter? Perhaps not: with a decent spin-up, the physics of the model should ensure the atmosphere comes into balance with its physical surroundings. But what if we are continously nudging the model with GFS data? Then it could be a problem.
The Web Mercartor Projection
Interestingly, the “Web Mercator” projection widely used in web mapping uses a similar scheme. The Earth is assumed to be spherical, yet the coordinates of points are required to be relative to the WGS84 ellipsoid. Wikipedia claims that this makes it slighly non-conformal, so presumably the same applies to LCC defined in WRF. How much affect this non-conformity has on the underlying equations in the model, I have no idea. I suspect it is insignificant compared to other sources of error.
Interpolating and reprojecting WRF output
There are two ways to interpolate WRF output to locations
- Use the latitude and longitude arrays as coordinates, and interpolate using them (ala rcm2points function in NCL)
- use the projection information to convert a latitide and longitiude into a grid-point index (ideally a decimal index), and then interpolate from the surrounding grid points (ala wrf_ll_to_ij in ncl)
In both of those situations, we shoudl get the correct result if we just take our geodetic WGS84 latitude and longitude and plug it in. However, if we want to reproject a whole WRF domain for the purpose of displaying it overlaid with other datasets, then we need to be a bit more cautious. There is a blog post here which deals with some of the same issues.
After a long winded search, I’ve come to conclusion that everything it probably OK. Observations are geolocated with a latitude relevant to the WGS94 ellipsoid, spherical models treat this as a spherical latitude. Most users do the same, and it probably all works.