Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

New functions for calculation of geopotential from height and vice versa #678

Merged
merged 9 commits into from Dec 18, 2017
Merged

New functions for calculation of geopotential from height and vice versa #678

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
ad7c444
Calculation of geopotential and geopot. height
Dec 8, 2017
696bbcb
New function for conversion from geopot to height
Dec 11, 2017
145ed1e
Update and style changes to geopotential functions
Dec 12, 2017
6e81070
Add tests for geopotential functions and fix typo
Dec 12, 2017
c9b3df9
Replace single with double quotes in docstrings
Dec 12, 2017
1d5ae99
Fix alphabetical import order
Dec 13, 2017
237e8e6
Fix lint error
Dec 13, 2017
6f90805
Fix the docstring test failure
Dec 13, 2017
51a9a6d
Parenthesis on line continuation in "constants"
Dec 18, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
73 changes: 72 additions & 1 deletion metpy/calc/basic.py
View file
Open in desktop
Expand Up @@ -16,7 +16,7 @@

import numpy as np

from ..constants import g, omega, Rd
from ..constants import g, omega, Rd, Re, G, me
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The checker doesn't like the order on these; it's alphabetization says:
"G, g, me, omega, Rd, Re"

from ..package_tools import Exporter
from ..units import atleast_1d, check_units, masked_array, units

Expand Down Expand Up @@ -257,6 +257,77 @@ def pressure_to_height_std(pressure):
return (t0 / gamma) * (1 - (pressure / p0).to('dimensionless')**(Rd * gamma / g))


@exporter.export
@check_units('[length]')
def height_to_geopotential(height):
r"""Computes geopotential for a given height in meters
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One of the typical styling things is to have this first line as an imperative and end in a period. Also, single quotes everywhere is what we've standardized on. Also, since we have the unit framework, we try to leave units out of the docstrings unless absolutely necessary. So can you change to:

'''Compute geopotential for a given height.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In basic.py all the docstrings contained double quotes so I replaced them in all the docstrings. Also changed to imperative and a period in the end.


Parameters
----------
height : `pint.Quantity`
Height above sea level (array_like)

Returns
-------
`pint.Quantity`
The corresponding geopotential value(s)

Examples
--------
>>> from metpy.constants import g,Re,G,me
>>> import metpy.calc
>>> from metpy.units import units
>>> height = np.linspace(0,10000, num = 11) * units.m
>>> geopot = metpy.calc.height_to_geopotential(height)
>>> geopot
<Quantity([ 0. 9817.70342881 19632.32592389 29443.86893527
39252.33391207 49057.72230251 58860.0355539 68659.27511263
78455.4424242 88248.53893319 98038.56608327],
'meter ** 2 / second ** 2')>
"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should be a blank line at the end of the docstring (inside the quotes).

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added the blank lines in the docstring(s).

# Calculate geopotential
geopot = G * me * ((1 / Re) - (1 / (Re + height)))

return geopot


@exporter.export
def geopotential_to_height(geopot):
r"""Computes height from a given geopotential
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here.


Parameters
----------
geopotential : `pint.Quantity`
Geopotential (array_like)

Returns
-------
`pint.Quantity`
The corresponding height value(s)

Examples
--------
>>> from metpy.constants import g,Re,G,me
>>> import metpy.calc
>>> from metpy.units import units
>>> height = np.linspace(0,10000, num = 11) * units.m
>>> geopot = metpy.calc.height_to_geopotential(height)
>>> geopot
<Quantity([ 0. 9817.70342881 19632.32592389 29443.86893527
39252.33391207 49057.72230251 58860.0355539 68659.27511263
78455.4424242 88248.53893319 98038.56608327],
'meter ** 2 / second ** 2')>
>>> height = metpy.calc.geopotential_to_height(geopot)
>>> height
<Quantity([ 0. 1000. 2000. 3000. 4000. 5000. 6000. 7000.
8000. 9000.10000.], 'meter')>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need some more separation spaces, especially between 9000. and 10000. So should be:

8000.  9000.  10000.

"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Blank line.

# Calculate geopotential
height = (((1 / Re) - (geopot / (G * me))) ** -1) - Re

return height


@exporter.export
@check_units('[length]')
def height_to_pressure_std(height):
Expand Down
26 changes: 15 additions & 11 deletions metpy/constants.py
View file
Open in desktop
Expand Up @@ -6,17 +6,19 @@

Earth
-----
======================== =============== =========== ========================== =======================================
Name Symbol Short Name Units Description
------------------------ --------------- ----------- -------------------------- ---------------------------------------
earth_avg_radius :math:`R_e` Re :math:`\text{m}` Avg. radius of the Earth
earth_gravity :math:`g` g :math:`\text{m s}^{-2}` Avg. gravity acceleration on Earth
earth_avg_angular_vel :math:`\Omega` omega :math:`\text{rad s}^{-1}` Avg. angular velocity of Earth
earth_sfc_avg_dist_sun :math:`d` d :math:`\text{m}` Avg. distance of the Earth from the Sun
earth_solar_irradiance :math:`S` S :math:`\text{W m}^{-2}` Avg. solar irradiance of Earth
earth_max_declination :math:`\delta` delta :math:`\text{degrees}` Max. solar declination angle of Earth
earth_orbit_eccentricity :math:`e` :math:`\text{None}` Avg. eccentricity of Earth's orbit
======================== =============== =========== ========================== =======================================
======================== =============== =========== ======================================= ======================================
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you need one more "=" at the end of the line.

Name Symbol Short Name Units Description
------------------------ --------------- ----------- --------------------------------------- --------------------------------------
earth_avg_radius :math:`R_e` Re :math:`\text{m}` Avg. radius of the Earth
earth_gravity :math:`g` g :math:`\text{m s}^{-2}` Avg. gravity acceleration on Earth
gravitational_constant :math:`G` G :math:`\text{m}^{3} {kg}^{-1} {s}^{-2}` Gravitational constant
earth_avg_angular_vel :math:`\Omega` omega :math:`\text{rad s}^{-1}` Avg. angular velocity of Earth
earth_sfc_avg_dist_sun :math:`d` d :math:`\text{m}` Avg. distance of the Earth from the Sun
earth_solar_irradiance :math:`S` S :math:`\text{W m}^{-2}` Avg. solar irradiance of Earth
earth_max_declination :math:`\delta` delta :math:`\text{degrees}` Max. solar declination angle of Earth
earth_orbit_eccentricity :math:`e` :math:`\text{None}` Avg. eccentricity of Earth's orbit
earth_mass :math:`m_e` me :math:`\text{kg}` Total mass of the Earth (approx)
======================== =============== =========== ======================================= =======================================

Water
-----
Expand Down Expand Up @@ -70,11 +72,13 @@
earth_gravity = g = units.Quantity(1.0, units.gravity).to('m / s^2')
# Taken from GEMPAK constants
Re = earth_avg_radius = 6.3712e6 * units.m
G = gravitational_constant = units.Quantity(6.674e-11, 'm^3 / kg / s^2')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you replace the hand-coded value here with the value from pint:

G = gravitational_constant = units.newtonian_constant_of_gravitation

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure if I've done this correctly, would be good to check!

omega = earth_avg_angular_vel = 2 * units.pi / units.sidereal_day
d = earth_sfc_avg_dist_sun = 1.496e11 * units.m
S = earth_solar_irradiance = units.Quantity(1.368e3, 'W / m^2')
delta = earth_max_declination = 23.45 * units.deg
earth_orbit_eccentricity = 0.0167
earth_mass = me = 5.9722e24 * units.kg

# molar gas constant
R = units.Quantity(1.0, units.R).to('J / K / mol')
Expand Down