geopulse.geo

Geographic projections used for line integrals of the geoelectric field.

Two projections are provided:

Both return (x_east_m, y_north_m) offsets in metres from a projection origin (lat0_deg, lon0_deg).

The WGS84 form is the default for PowerGridNetwork and PipelineNetwork so that GIC benchmark results match published papers to sub-percent precision.

Notes

Horton (2012) Appendix — the north-south and east-west local distances between two lat/lon points A, B in kilometres are:

\[\begin{split}\phi &= (\mathrm{lat}_A + \mathrm{lat}_B) / 2 \\ L_N &= (111.133 - 0.56\,\cos 2\phi)\,\Delta\mathrm{lat} \\ L_E &= (111.5065 - 0.1872\,\cos 2\phi)\,\cos\phi\,\Delta\mathrm{lon}\end{split}\]

with Δlat, Δlon in degrees. This is the compact form of the exact meridian- and prime-vertical-radius-of-curvature expressions:

\[\begin{split}M(\phi) &= \frac{a(1 - e^2)}{(1 - e^2 \sin^2\phi)^{3/2}} \\ N(\phi) &= \frac{a}{\sqrt{1 - e^2 \sin^2\phi}} \\ L_N &= \frac{\pi}{180}\,M(\phi)\,\Delta\mathrm{lat} \\ L_E &= \frac{\pi}{180}\,N(\phi)\,\cos\phi\,\Delta\mathrm{lon}\end{split}\]

The implementation uses the exact M, N form so we stay accurate at polar latitudes where the truncated series in the paper’s Appendix starts to lose digits.

References

Functions

latlon_to_local_xy_spherical_m(lat_deg, ...)

Equirectangular projection about (lat0, lon0) — spherical Earth.

latlon_to_local_xy_wgs84_m(lat_deg, lon_deg, ...)

Equirectangular projection about (lat0, lon0) — WGS84 ellipsoid.

meridian_radius_m(lat_deg)

WGS84 meridional radius of curvature M(φ) in metres.

prime_vertical_radius_m(lat_deg)

WGS84 prime-vertical radius of curvature N(φ) in metres.
geopulse.geo.latlon_to_local_xy_spherical_m(lat_deg, lon_deg, lat0_deg, lon0_deg)[source]

Equirectangular projection about (lat0, lon0) — spherical Earth.

Uses R_EARTH_M (mean radius 6371 km). Cheap and adequate to ~0.3 % for a few hundred kilometres, but drifts at continent scale. Prefer latlon_to_local_xy_wgs84_m() for anything published.

Parameters:

  • lat_deg (float) – Point coordinates in degrees.

  • lon_deg (float) – Point coordinates in degrees.

  • lat0_deg (float) – Projection origin in degrees.

  • lon0_deg (float) – Projection origin in degrees.

Return type:

tuple[float, float]

Returns:

x_east_m, y_north_m (float) – Eastward and northward offsets in metres.

geopulse.geo.latlon_to_local_xy_wgs84_m(lat_deg, lon_deg, lat0_deg, lon0_deg)[source]

Equirectangular projection about (lat0, lon0) — WGS84 ellipsoid.

Uses meridian- and prime-vertical-radii-of-curvature at the mean latitude φ = (lat + lat0) / 2 to convert degree offsets into east/north offsets in metres. Matches Horton (2012) Appendix A3/A7 exactly (the 111.133 − 0.56 cos 2φ etc. truncated forms are the small-e² series expansion of this same expression).

Parameters:

  • lat_deg (float) – Point coordinates in degrees.

  • lon_deg (float) – Point coordinates in degrees.

  • lat0_deg (float) – Projection origin in degrees.

  • lon0_deg (float) – Projection origin in degrees.

Return type:

tuple[float, float]

Returns:

x_east_m, y_north_m (float) – Eastward and northward offsets in metres.

Examples

>>> # 1° latitude at 45° N ≈ 111.13 km per WGS84.
>>> x, y = latlon_to_local_xy_wgs84_m(46.0, -75.0, 45.0, -75.0)
>>> 111_000 < y < 112_000
True
geopulse.geo.meridian_radius_m(lat_deg)[source]

WGS84 meridional radius of curvature M(φ) in metres.

M(φ) = a(1 − e²) / (1 − e² sin²φ)^{3/2} — the radius of curvature of the ellipsoid along a meridian (north-south).

Return type:

float

Parameters:

lat_deg (float)

geopulse.geo.prime_vertical_radius_m(lat_deg)[source]

WGS84 prime-vertical radius of curvature N(φ) in metres.

N(φ) = a / √(1 − e² sin²φ) — the radius of curvature in the plane perpendicular to a meridian, at latitude φ.

Return type:

float

Parameters:

lat_deg (float)