pycontrails.physics.geo¶

Tools for spherical geometry, solar radiation, and wind advection.

Functions

`advect_horizontal`(longitude, latitude, ...)	Advect a particle in the horizontal plane.
`advect_latitude`(latitude, v_wind, dt)	Calculate the latitude of a particle after time `dt` caused by advection due to wind.
`advect_level`(level, vertical_velocity, ...)	Calculate the pressure level of a particle after time `dt`.
`advect_longitude`(longitude, latitude, u_wind, dt)	Calculate the longitude of a particle after time dt caused by advection due to wind.
`advect_longitude_and_latitude_near_poles`(...)	Advect a particle near the poles.
`azimuth`(lons0, lats0, lons1, lats1)	Calculate angle relative to true north for set of coordinates.
`azimuth_to_direction`(azimuth_, latitude)	Calculate rectangular direction from spherical azimuth.
`cosine_solar_zenith_angle`(longitude, ...)	Calculate the cosine of the solar zenith angle.
`days_since_reference_year`(time[, ref_year])	Calculate the days elapsed since the start of the reference year.
`domain_surface_area`([spatial_bbox, ...])	Calculate surface area in the provided spatial bounding box.
`forward_azimuth`(lons, lats, az, dist)	Calculate coordinates along forward azimuth.
`grid_surface_area`(longitude, latitude)	Calculate surface area that is covered by each pixel in a longitude-latitude grid.
`haversine`(lons0, lats0, lons1, lats1)	Calculate haversine distance between points in (lons0, lats0) and (lons1, lats1).
`hours_since_start_of_day`(time)	Calculate the hours elapsed since the start of day (00:00:00 UTC).
`longitudinal_angle`(lons0, lats0, lons1, lats1)	Calculate angle with longitudinal axis for sequence of segments.
`orbital_correction_for_solar_hour_angle`(...)	Calculate correction to the solar hour angle due to Earth's orbital location.
`orbital_position`(time)	Calculate the orbital position of Earth to a reference point set at the start of year.
`segment_angle`(longitude, latitude)	Calculate the angle between coordinate segments and the longitudinal axis.
`segment_azimuth`(longitude, latitude)	Calculate the angle between coordinate segments and true north.
`segment_haversine`(longitude, latitude)	Calculate haversine distance between consecutive points along path.
`segment_length`(longitude, latitude, altitude)	Calculate the segment length between coordinates by assuming a great circle distance.
`solar_constant`(theta_rad)	Calculate the solar electromagnetic radiation per unit area from orbital position.
`solar_declination_angle`(theta_rad)	Calculate the solar declination angle from the orbital position in radians (theta_rad).
`solar_direct_radiation`(longitude, latitude, time)	Calculate the instantaneous theoretical solar direct radiation (SDR).
`solar_hour_angle`(longitude, time, theta_rad)	Calculate the sun's East to West angular displacement around the polar axis.
`spatial_bounding_box`(longitude, latitude[, ...])	Construct rectangular spatial bounding box from a set of waypoints.

pycontrails.physics.geo.advect_horizontal(longitude, latitude, u_wind, v_wind, dt)¶

Advect a particle in the horizontal plane.

This function calls advect_longitude() and advect_latitude() when the position is far from the poles (<= 80.0 degrees). When the position is near the poles (> 80.0 degrees), advect_longitude_and_latitude_near_poles() is used instead.

Parameters:

longitude (npt.NDArray[np.floating]) – Original longitude, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – Original latitude, [\(\deg\)]
u_wind (npt.NDArray[np.floating]) – Wind speed in the longitudinal direction, [\(m s^{-1}\)]
v_wind (npt.NDArray[np.floating]) – Wind speed in the latitudinal direction, [\(m s^{-1}\)]
dt (npt.NDArray[np.timedelta64] | np.timedelta64) – Advection timestep

Returns:

tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]] – New longitude and latitude values, [\(\deg\)]

pycontrails.physics.geo.advect_latitude(latitude, v_wind, dt)¶

Calculate the latitude of a particle after time dt caused by advection due to wind.

Note

It is possible for advected latitude values to lie outside of the WGS84 domain [-90, 90]. In Cocip models, latitude values close to the poles create an end of life condition, thereby avoiding this issue. In practice, such situations are very rare.

These polar divergence issues could also be addressed by reflecting the longitude values 180 degrees via a spherical equivalence such as (lon, lat) ~ (lon + 180, 180 - lat). This approach is not currently taken.

Parameters:

latitude (ArrayLike) – Original latitude, [\(\deg\)]
v_wind (ArrayLike) – Wind speed in the latitudinal direction, [\(m s^{-1}\)]
dt (numpy.ndarray) – Advection time delta

Returns:

ArrayLike – New latitude value, [\(\deg\)]

pycontrails.physics.geo.advect_level(level, vertical_velocity, rho_air, terminal_fall_speed, dt)¶

Calculate the pressure level of a particle after time dt.

This function calculates the new pressure level of a particle as a result of vertical advection caused by the vertical velocity and terminal fall speed.

Parameters:

level (ArrayLike) – Pressure level, [\(hPa\)]
vertical_velocity (ArrayLike) – Vertical velocity, [\(Pa s^{-1}\)]
rho_air (ArrayLike | float) – Air density, [\(kg m^{-3}\)]
terminal_fall_speed (ArrayLike | float) – Terminal fall speed of the particle, [\(m s^{-1}\)]
dt (npt.NDArray[np.timedelta64] | np.timedelta64) – Time delta for each waypoint

Returns:

ArrayLike – New pressure level, [\(hPa\)]

pycontrails.physics.geo.advect_longitude(longitude, latitude, u_wind, dt)¶

Calculate the longitude of a particle after time dt caused by advection due to wind.

Automatically wrap over the antimeridian if necessary.

Parameters:

longitude (ArrayLike) – Original longitude, [\(\deg\)]
latitude (ArrayLike) – Original latitude, [\(\deg\)]
u_wind (ArrayLike) – Wind speed in the longitudinal direction, [\(m s^{-1}\)]
dt (numpy.ndarray) – Advection timestep

Returns:

ArrayLike – New longitude value, [\(\deg\)]

pycontrails.physics.geo.advect_longitude_and_latitude_near_poles(longitude, latitude, u_wind, v_wind, dt)¶

Advect a particle near the poles.

This function calculates the longitude and latitude of a particle after time dt caused by advection due to wind near the poles (above 80 degrees North and South).

Automatically wrap over the antimeridian if necessary.

Parameters:

longitude (npt.NDArray[np.floating]) – Original longitude, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – Original latitude, [\(\deg\)]
u_wind (npt.NDArray[np.floating]) – Wind speed in the longitudinal direction, [\(m s^{-1}\)]
v_wind (npt.NDArray[np.floating]) – Wind speed in the latitudinal direction, [\(m s^{-1}\)]
dt (npt.NDArray[np.timedelta64] | np.timedelta64) – Advection timestep

Returns:

tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]] – New longitude and latitude values, [\(\deg\)]

Notes

Near the poles, the longitude and latitude is converted to a 2-D Cartesian-like coordinate system to avoid numerical instabilities and singularities caused by convergence of meridians.

pycontrails.physics.geo.azimuth(lons0, lats0, lons1, lats1)¶

Calculate angle relative to true north for set of coordinates.

Parameters:

lons0 (npt.NDArray[np.floating]) – Longitude values of initial endpoints, [\(\deg\)].
lats0 (npt.NDArray[np.floating]) – Latitude values of initial endpoints, [\(\deg\)].
lons1 (npt.NDArray[np.floating]) – Longitude values of terminal endpoints, [\(\deg\)].
lats1 (npt.NDArray[np.floating]) – Latitude values of terminal endpoints, [\(\deg\)].

References

[Wikipedia contributors, 2023]

Returns:: npt.NDArray[np.floating] – Azimuth relative to true north (\(0\deg\)), [\(\deg\)]

See also

longitudinal_angle()

pycontrails.physics.geo.azimuth_to_direction(azimuth_, latitude)¶

Calculate rectangular direction from spherical azimuth.

This implementation uses the equation

cos(latitude) / tan(azimuth) = sin_a / cos_a

to solve for sin_a and cos_a.

Parameters:

azimuth_ (npt.NDArray[np.floating]) – Angle measured clockwise from true north, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – Latitude value of the point, [\(\deg\)]

Returns:

tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]] – A tuple of sine and cosine values.

pycontrails.physics.geo.cosine_solar_zenith_angle(longitude, latitude, time, theta_rad)¶

Calculate the cosine of the solar zenith angle.

Return (\(\cos(\theta)\)), where \(\theta\) is the angle between the sun and the vertical direction.

Parameters:

longitude (ArrayLike) – Longitude, [\(\deg\)]
latitude (ArrayLike) – Latitude, [\(\deg\)]
time (ArrayLike) – Time, formatted as np.datetime64
theta_rad (ArrayLike) – Orbital position, [\(rad\)]. Output of orbital_position().

Returns:

ArrayLike – Cosine of the solar zenith angle

References

[Wikipedia contributors, 2023]

pycontrails.physics.geo.days_since_reference_year(time, ref_year=2000)¶

Calculate the days elapsed since the start of the reference year.

Parameters:

time (ArrayLike) – ArrayLike of np.datetime64 times
ref_year (int, optional) – Year of reference

Returns:

ArrayLike – Days elapsed since the reference year. Output dtype is np.float64.

Raises:

RuntimeError – Raises when reference year is greater than the time of time element

pycontrails.physics.geo.domain_surface_area(spatial_bbox=(-180.0, -90.0, 180.0, 90.0), spatial_grid_res=0.5)¶

Calculate surface area in the provided spatial bounding box.

Parameters:

spatial_bbox (tuple[float, float, float, float]) – Spatial bounding box, (lon_min, lat_min, lon_max, lat_max), [\(\deg\)]
spatial_grid_res (float) – Spatial grid resolution, [\(\deg\)]

Returns:

float – Domain surface area, [\(m^{2}\)]

pycontrails.physics.geo.forward_azimuth(lons, lats, az, dist)¶

Calculate coordinates along forward azimuth.

This function is identical to the pyproj.Geod.fwd method when working on a spherical earth. Both signatures are also identical. This implementation is generally more performant.

Parameters:

lons (npt.NDArray[np.floating]) – Array of longitude values.
lats (npt.NDArray[np.floating]) – Array of latitude values.
az (npt.NDArray[np.floating] | float) – Azimuth, measured in [\(\deg\)].
dist (npt.NDArray[np.floating] | float) – Distance [\(m\)] between initial longitude latitude values and point to be computed.

Returns:

tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]] – Tuple of longitude latitude arrays.

See also

:meth:pyproj.Geod.fwd

pycontrails.physics.geo.grid_surface_area(longitude, latitude)¶

Calculate surface area that is covered by each pixel in a longitude-latitude grid.

Parameters:

longitude (npt.NDArray[np.floating]) – Longitude coordinates in a longitude-latitude grid, [\(\deg\)]. Must be in ascending order.
latitude (npt.NDArray[np.floating]) – Latitude coordinates in a longitude-latitude grid, [\(\deg\)]. Must be in ascending order.

Returns:

xarray.DataArray – Surface area of each pixel in a longitude-latitude grid, [\(m^{2}\)]

References

https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2004/msg00023.html

pycontrails.physics.geo.haversine(lons0, lats0, lons1, lats1)¶

Calculate haversine distance between points in (lons0, lats0) and (lons1, lats1).

Handles coordinates crossing the antimeridian line (-180, 180).

Parameters:

lons0, lats0 (ArrayLike) – Coordinates of initial points, [\(\deg\)]
lons1, lats1 (ArrayLike) – Coordinates of terminal points, [\(\deg\)]

Returns:

ArrayLike – Distances between corresponding points. [\(m\)]

Notes

This formula does not take into account the non-spheroidal (ellipsoidal) shape of the Earth. Originally referenced from https://andrew.hedges.name/experiments/haversine/.

References

[Calculate Distance and Bearing between Two Latitude/Longitude Points Using Haversine Formula in JavaScript, n.d.]

See also

sklearn.metrics.pairwise.haversine_distances(): Compute the Haversine distance
pyproj.Geod: Performs forward and inverse geodetic, or Great Circle, computations

pycontrails.physics.geo.hours_since_start_of_day(time)¶

Calculate the hours elapsed since the start of day (00:00:00 UTC).

Parameters:: time (ArrayLike) – ArrayLike of np.datetime64 times
Returns:: ArrayLike – Hours elapsed since the start of today day. Output dtype is np.float64.

pycontrails.physics.geo.longitudinal_angle(lons0, lats0, lons1, lats1)¶

Calculate angle with longitudinal axis for sequence of segments.

Parameters:

lons0 (npt.NDArray[np.floating]) – Longitude values of initial endpoints, [\(\deg\)].
lats0 (npt.NDArray[np.floating]) – Latitude values of initial endpoints, [\(\deg\)].
lons1 (npt.NDArray[np.floating]) – Longitude values of terminal endpoints, [\(\deg\)].
lats1 (npt.NDArray[np.floating]) – Latitude values of terminal endpoints, [\(\deg\)].

References

[Wikipedia contributors, 2023]

Returns:

sin_a (npt.NDArray[np.floating]) – Sine values.
cos_a (npt.NDArray[np.floating]) – Cosine values.

pycontrails.physics.geo.orbital_correction_for_solar_hour_angle(theta_rad)¶

Calculate correction to the solar hour angle due to Earth’s orbital location.

Parameters:: theta_rad (ArrayLike) – Orbital position, [\(rad\)]
Returns:: ArrayLike – Correction to the solar hour angle as a result of Earth’s orbital location, [\(\deg\)]

References

[Paltridge et al., 1976]

Notes

Tested against [NOAA, n.d.]

pycontrails.physics.geo.orbital_position(time)¶

Calculate the orbital position of Earth to a reference point set at the start of year.

Parameters:: time (ArrayLike) – ArrayLike of np.datetime64 times
Returns:: ArrayLike – Orbital position of Earth, [\(rad\)]

pycontrails.physics.geo.segment_angle(longitude, latitude)¶

Calculate the angle between coordinate segments and the longitudinal axis.

np.nan is added to the final value so the length of the output is the same as the inputs.

Parameters:

longitude (npt.NDArray[np.floating]) – Longitude values, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – Latitude values, [\(\deg\)]

Returns:

tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]] – sin(a), cos(a), where a is the angle between the segment and the longitudinal axis. Final entry of each array is set to np.nan.

References

[Wikipedia contributors, 2023]

Notes

        (lon_2, lat_2)  X
                       /|
                      / |
                     /  |
                    /   |
                   /    |
                  /     |
                 /      |
(lon_1, lat_1)  X -------> longitude (x-axis)

See also

longitudinal_angle()

pycontrails.physics.geo.segment_azimuth(longitude, latitude)¶

Calculate the angle between coordinate segments and true north.

np.nan is added to the final value so the length of the output is the same as the inputs.

Parameters:

longitude (npt.NDArray[np.floating]) – Longitude values, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – Latitude values, [\(\deg\)]

Returns:

npt.NDArray[np.floating] – Azimuth relative to true north (\(0\deg\)), [\(\deg\)] Final entry of each array is set to np.nan.

References

[Wikipedia contributors, 2023]

See also

azimuth()

pycontrails.physics.geo.segment_haversine(longitude, latitude)¶

Calculate haversine distance between consecutive points along path.

Parameters:

longitude (npt.NDArray[np.floating]) – 1D Longitude values with index corresponding to latitude inputs, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – 1D Latitude values with index corresponding to longitude inputs, [\(\deg\)]

Returns:

npt.NDArray[np.floating] – Haversine distance between (lat_i, lon_i) and (lat_i+1, lon_i+1), [\(m\)] The final entry of the output is set to nan.

See also

pyproj.Geod.line_lengths()

pycontrails.physics.geo.segment_length(longitude, latitude, altitude)¶

Calculate the segment length between coordinates by assuming a great circle distance.

Requires coordinates to be in EPSG:4326. Lengths are calculated using both horizontal and vertical displacement of segments.

np.nan is added to the final value so the length of the output is the same as the inputs.

Parameters:

longitude (npt.NDArray[np.floating]) – Longitude values, [\(\deg\)]
latitude (npt.NDArray[np.floating]) – Latitude values, [\(\deg\)]
altitude (npt.NDArray[np.floating]) – Altitude values, [\(m\)]

Returns:

npt.NDArray[np.floating] – Array of distances in [\(m\)] between coordinates. Final entry of each array is set to np.nan.

See also

haversine(), segment_haversine()

pycontrails.physics.geo.solar_constant(theta_rad)¶

Calculate the solar electromagnetic radiation per unit area from orbital position.

On average, the extraterrestrial irradiance is 1367 W/m**2 and varies by +- 3% as the Earth orbits the sun.

Parameters:: theta_rad (ArrayLike) – Orbital position, [\(rad\)]. Use orbital_position() to calculate the orbital position from time input.
Returns:: ArrayLike – Solar constant, [\(W m^{-2}\)]

References

Notes

\(orbital_effect = (R_{av} / R)^{2}\) where \(R\) is the separation of Earth from the sun and \(R_{av}\) is the mean separation.

pycontrails.physics.geo.solar_declination_angle(theta_rad)¶

Calculate the solar declination angle from the orbital position in radians (theta_rad).

The solar declination angle is the angle between the rays of the Sun and the plane of the Earth’s equator.

It has a range of between -23.5 (winter solstice) and +23.5 (summer solstice) degrees.

Parameters:: theta_rad (ArrayLike) – Orbital position, [\(rad\)]. Output of orbital_position().
Returns:: ArrayLike – Solar declination angle, [\(\deg\)]

References

[Paltridge et al., 1976]

Notes

Tested against [NOAA, n.d.]

pycontrails.physics.geo.solar_direct_radiation(longitude, latitude, time, threshold_cos_sza=0.0)¶

Calculate the instantaneous theoretical solar direct radiation (SDR).

Parameters:

longitude (ArrayLike) – Longitude, [\(\deg\)]
latitude (ArrayLike) – Latitude, [\(\deg\)]
time (ArrayLike) – Time, formatted as np.datetime64
threshold_cos_sza (float, optional) – Set the SDR to 0 when the cosine_solar_zenith_angle() is below a certain value. By default, set to 0.

Returns:

ArrayLike – Solar direct radiation of incoming radiation, [\(W m^{-2}\)]

References

[Laboratory, 2022]

pycontrails.physics.geo.solar_hour_angle(longitude, time, theta_rad)¶

Calculate the sun’s East to West angular displacement around the polar axis.

The solar hour angle is an expression of time in angular measurements: the value of the hour angle is zero at noon, negative in the morning, and positive in the afternoon, increasing by 15 degrees per hour.

Parameters:

longitude (ArrayLike) – Longitude, [\(\deg\)]
time (ArrayLike) – ArrayLike of np.datetime64 times
theta_rad (ArrayLike) – Orbital position, [\(rad\)]. Output of orbital_position().

Returns:

ArrayLike – Solar hour angle, [\(\deg\)]

pycontrails.physics.geo.spatial_bounding_box(longitude, latitude, buffer=1.0)¶

Construct rectangular spatial bounding box from a set of waypoints.

Parameters:

longitude (numpy.ndarray) – 1D Longitude values with index corresponding to longitude inputs, [\(\deg\)]
latitude (numpy.ndarray) – 1D Latitude values with index corresponding to latitude inputs, [\(\deg\)]
buffer (float) – Add buffer to rectangular spatial bounding box, [\(\deg\)]

Returns:

tuple[float, float, float, float] – Spatial bounding box, (lon_min, lat_min, lon_max, lat_max), [\(\deg\)]

Examples

>>> rng = np.random.default_rng(654321)
>>> lon = rng.uniform(-180, 180, size=30)
>>> lat = rng.uniform(-90, 90, size=30)
>>> spatial_bounding_box(lon, lat)
(np.float64(-168.0), np.float64(-77.0), np.float64(155.0), np.float64(82.0))