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

Implement gwfs function from GEMPAK #911

Merged
merged 12 commits into from
Aug 6, 2018
Merged

Implement gwfs function from GEMPAK #911

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
77782e7
Fix non-tuple NumPy indexing in MetPy
jthielen Aug 1, 2018
764b722
Add GWFS function from GEMPAK
sgdecker Aug 2, 2018
841c8ea
Add tests and docs for gwfs()
sgdecker Aug 2, 2018
b1d7f0e
Fix whitespace issues re pep8
sgdecker Aug 2, 2018
d531f34
Merge pull request #910 from jthielen/non-tuple-indexing-cleanup
dopplershift Aug 2, 2018
261188f
Add GWFS function from GEMPAK
sgdecker Aug 2, 2018
10418d7
Add tests and docs for gwfs()
sgdecker Aug 2, 2018
027bb97
Fix whitespace issues re pep8
sgdecker Aug 2, 2018
4f898e2
Extend functionality of gwfs and update docs
sgdecker Aug 3, 2018
067d33b
Rename gwfs() to smooth_gaussian()
sgdecker Aug 3, 2018
582237d
Merge branch 'add_gwf_funs' of https://github.com/sgdecker/metpy into…
sgdecker Aug 3, 2018
2fdb0d4
Update GEMPAK Conversion Guide
sgdecker Aug 3, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions docs/gempak.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -469,11 +469,11 @@ blue is uncertain of parity, and white is unevaluated.
<td></td>
</tr>
<tr>
<td>GWFS(S, N)</td>
<td>Filter with normal distribution of weights</td>
<td></td>
<td></td>
<td></td>
<td class="tg-implemented">GWFS(S, N)</td>
<td class="tg-implemented">Filter with normal distribution of weights</td>
<td class="tg-implemented"><a href="api/generated/metpy.calc.gwfs.html#metpy.calc.gwfs">metpy.calc.gwfs</a></td>
<td class="tg-yes">Yes</td>
<td class="tg-yes">Yes</td>
<td></td>
</tr>
<tr>
Expand Down
85 changes: 85 additions & 0 deletions metpy/calc/basic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,8 @@

import numpy as np

from scipy.ndimage import gaussian_filter

from ..constants import G, g, me, omega, Rd, Re
from ..deprecation import deprecated
from ..package_tools import Exporter
Expand Down Expand Up @@ -614,6 +616,89 @@ def sigma_to_pressure(sigma, psfc, ptop):
return sigma * (psfc - ptop) + ptop


@exporter.export
@preprocess_xarray
def gwfs(scalar_grid, n):
Copy link
Member

Choose a reason for hiding this comment

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

We need a name for this--replicating GEMPAK functionality is a big 👍 . Replicating naming from the 80's when you had to pay per-character (or so it seems) is something I'd like to not replicate. 😁

How about gaussian_smooth?

Also, is n the most sensible argument for this function? Personally, I'd like to be able to specify the kernel width in some kind of physical units as well, but we don't have to do this right now. If this is the most sensible interface to you the user, then I'm happy to just roll with it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, gwfs is more of a placeholder name to make it clear what GEMPAK function it corresponds to. Looking ahead to eventual implementations of GEMPAK's sm5s and sm9s, maybe all of these should start with 'smooth', so smooth_gaussian, smooth_5pt, etc.?

Regarding n, it sounds like you are wanting the related GEMPAK function RDFS (Resolution Dependent Filter for Scalar), which takes dx instead of n as the second argument, and then behind the scenes computes an appropriate n based on dx and calls GWFS. But I suppose just because they are two separate GEMPAK functions doesn't mean they have to be separate MetPy functions.

Copy link
Member

Choose a reason for hiding this comment

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

smooth_gaussian sounds good to me.

If both GWFS and RDFS are doing the same calculation, just tweaking the parameter, I’d lean towards putting all of the functionality in smooth_gaussian, and just look at the arguments passed in to determine the proper course of action. That doesn’t have to be implemented in this PR, though—I’d be happy just to merge the GWFS equivalent.

"""Filter with normal distribution of weights.

Parameters
----------
scalar_grid : `pint.Quantity`
Some two-dimensional scalar grid

n : int
Degree of filtering

Returns
-------
`pint.Quantity`
The filtered 2D scalar grid

Copy link

Choose a reason for hiding this comment

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

W293 blank line contains whitespace

Notes
-----
This function is a close replication of the GEMPAK function GWFS,
but is not identical. The following notes are incorporated from
the GEMPAK source code:

This function smoothes a scalar grid using a moving average
low-pass filter whose weights are determined by the normal
(Gaussian) probability distribution function for two dimensions.
The weight given to any grid point within the area covered by the
moving average for a target grid point is proportional to

EXP [ -( D ** 2 ) ],

where D is the distance from that point to the target point divided
by the standard deviation of the normal distribution. The value of
the standard deviation is determined by the degree of filtering
requested. The degree of filtering is specified by an integer.
This integer is the number of grid increments from crest to crest
of the wave for which the theoretical response is 1/e = .3679. If
the grid increment is called delta_x, and the value of this integer
is represented by N, then the theoretical filter response function
value for the N * delta_x wave will be 1/e. The actual response
function will be greater than the theoretical value.

The larger N is, the more severe the filtering will be, because the
response function for all wavelengths shorter than N * delta_x
will be less than 1/e. Furthermore, as N is increased, the slope
of the filter response function becomes more shallow; so, the
response at all wavelengths decreases, but the amount of decrease
lessens with increasing wavelength. (The theoretical response
function can be obtained easily--it is the Fourier transform of the
weight function described above.)

The area of the patch covered by the moving average varies with N.
As N gets bigger, the smoothing gets stronger, and weight values
farther from the target grid point are larger because the standard
deviation of the normal distribution is bigger. Thus, increasing
N has the effect of expanding the moving average window as well as
changing the values of weights. The patch is a square covering all
points whose weight values are within two standard deviations of the
mean of the two dimensional normal distribution.

The key difference between GEMPAK's GWFS and this function is that,
in GEMPAK, the leftover weight values representing the fringe of the
distribution are applied to the target grid point. In this
function, the leftover weights are not used.

When this function is invoked, the first argument is the grid to be
smoothed, the second is the value of N as described above:

GWFS ( S, N )

where N > 1. If N <= 1, N = 2 is assumed. For example, if N = 4,
then the 4 delta x wave length is passed with approximate response
1/e.
Copy link

Choose a reason for hiding this comment

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

W291 trailing whitespace

"""
n = int(round(n))
if n < 2:
n = 2
sgma = n / (2*np.pi)
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

res = gaussian_filter(scalar_grid, sgma, truncate=2*np.sqrt(2))
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

return res


Copy link

Choose a reason for hiding this comment

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

W293 blank line contains whitespace

def _check_radians(value, max_radians=2 * np.pi):
"""Input validation of values that could be in degrees instead of radians.

Expand Down
47 changes: 46 additions & 1 deletion metpy/calc/tests/test_basic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

from metpy.calc import (add_height_to_pressure, add_pressure_to_height, apparent_temperature,
coriolis_parameter, geopotential_to_height, get_wind_components,
get_wind_dir, get_wind_speed, heat_index, height_to_geopotential,
get_wind_dir, get_wind_speed, gwfs, heat_index, height_to_geopotential,
height_to_pressure_std, pressure_to_height_std, sigma_to_pressure,
wind_components, wind_direction, wind_speed, windchill)
from metpy.deprecation import MetpyDeprecationWarning
Expand Down Expand Up @@ -436,3 +436,48 @@ def test_apparent_temperature_windchill():
truth = -18.9357 * units.degC
res = apparent_temperature(temperature, rel_humidity, wind)
assert_almost_equal(res, truth, 0)


def test_gwfs():
"""Test the gwfs function with a larger n."""
m = 10
s = np.zeros((m, m))
for i in np.ndindex(s.shape):
s[i] = i[0] + i[1]**2
s = gwfs(s, 4)
s_true = np.array([[ 0.40077472, 1.59215426, 4.59665817, 9.59665817, 16.59665817,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

25.59665817, 36.59665817, 49.59665817, 64.51108392, 77.87487258],
[ 1.20939518, 2.40077472, 5.40527863, 10.40527863, 17.40527863,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

26.40527863, 37.40527863, 50.40527863, 65.31970438, 78.68349304],
[ 2.20489127, 3.39627081, 6.40077472, 11.40077472, 18.40077472,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

27.40077472, 38.40077472, 51.40077472, 66.31520047, 79.67898913],
[ 3.20489127, 4.39627081, 7.40077472, 12.40077472, 19.40077472,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

28.40077472, 39.40077472, 52.40077472, 67.31520047, 80.67898913],
[ 4.20489127, 5.39627081, 8.40077472, 13.40077472, 20.40077472,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

29.40077472, 40.40077472, 53.40077472, 68.31520047, 81.67898913,],
Copy link

Choose a reason for hiding this comment

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

E231 missing whitespace after ','

[ 5.20489127, 6.39627081, 9.40077472, 14.40077472, 21.40077472,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

30.40077472, 41.40077472, 54.40077472, 69.31520047, 82.67898913],
[ 6.20489127, 7.39627081, 10.40077472, 15.40077472, 22.40077472,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

31.40077472, 42.40077472, 55.40077472, 70.31520047, 83.67898913],
[ 7.20489127, 8.39627081, 11.40077472, 16.40077472, 23.40077472,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

32.40077472, 43.40077472, 56.40077472, 71.31520047, 84.67898913],
[ 8.20038736, 9.3917669, 12.39627081, 17.39627081, 24.39627081,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['
E241 multiple spaces after ','

33.39627081, 44.39627081, 57.39627081, 72.31069656, 85.67448522],
[ 9.00900782, 10.20038736, 13.20489127, 18.20489127, 25.20489127,
Copy link

Choose a reason for hiding this comment

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

E201 whitespace after '['

34.20489127, 45.20489127, 58.20489127, 73.11931702, 86.48310568]])
assert_array_almost_equal(s, s_true)

def test_gwfs_small_n():
Copy link

Choose a reason for hiding this comment

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

E302 expected 2 blank lines, found 1

"""Test the gwfs function with a smaller n."""
m = 5
s = np.zeros((m, m))
for i in np.ndindex(s.shape):
s[i] = i[0] + i[1]**2
s = gwfs(s, 1)
s_true = [[0.0141798077, 1.02126971, 4.02126971, 9.02126971, 15.9574606],
[1.00708990, 2.01417981, 5.01417981, 10.0141798, 16.9503707],
[2.00708990, 3.01417981, 6.01417981, 11.0141798, 17.9503707],
[3.00708990, 4.01417981, 7.01417981, 12.0141798, 18.9503707],
[4.00000000, 5.00708990, 8.00708990, 13.0070899, 19.9432808]]
print(s-s_true)
Copy link

Choose a reason for hiding this comment

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

E226 missing whitespace around arithmetic operator

assert_array_almost_equal(s, s_true)