cross_section() - Code Metrics - Inspection of "extract_cross_section" - Unidata/MetPy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Pull Request — master (#518)

unknown

created 2017-07-25 17:22 UTC

cross_section() B

↳ Parent: Project

Complexity

Conditions

Size

Total Lines

104

Duplication

Lines	0
Ratio	0 %

Importance

Changes	1
Bugs	0	Features	0

Metric	Value
cc	2
c	1
b	0
f	0
dl	0
loc	104
rs	8.2857

How to fix Long Method

# Copyright (c) 2008-2017 MetPy Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
"""
Isentropic Cross-section
========================

Make an isentropic cross section with MetPy

With the function `metpy.calc.extract_cross_section`, data can be extracted from 3-dimensional
input along a specified line and plotted in a 2-dimensional cross-section.

.. image:: ../../_static/cross_section.gif

"""

################################
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import ipywidgets as widgets
from ipywidgets import interact_manual
import matplotlib.pyplot as plt
from netCDF4 import Dataset, num2date
import numpy as np

import metpy.calc as mcalc
from metpy.cbook import get_test_data
from metpy.units import units

#######################################
# **Getting the data**
#
# In this example, NARR reanalysis data for 18 UTC 04 April 1987 from the National Centers
# for Environmental Information (https://nomads.ncdc.noaa.gov) will be used.

data = Dataset(get_test_data('narr_example.nc', False))

##########################
print(list(data.variables))

#############################
# We will reduce the dimensionality of the data as it is pulled in to remove an empty time
# dimension. Additionally, units are required for input data, so the proper units will also
# be attached.


# Assign data to variable names
dtime = data.variables['Geopotential_height'].dimensions[0]
dlev = data.variables['Geopotential_height'].dimensions[1]
lat = data.variables['lat'][:]
lon = data.variables['lon'][:]
lev = data.variables[dlev][:] * units(data.variables[dlev].units)
times = data.variables[dtime]
vtimes = num2date(times[:], times.units)

temps = data.variables['Temperature']
tmp = temps[0, :] * units.kelvin
uwnd = data.variables['u_wind'][0, :] * units(data.variables['u_wind'].units)
vwnd = data.variables['v_wind'][0, :] * units(data.variables['v_wind'].units)
hgt = data.variables['Geopotential_height'][0, :] * units.meter
spech = (data.variables['Specific_humidity'][0, :] *
         units(data.variables['Specific_humidity'].units))


#########################
# Create the coordinate system and projection for plotting the inset map
crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)
tlatlons = crs.transform_points(ccrs.PlateCarree(), lon, lat)
tlons = tlatlons[:, :, 0]
tlats = tlatlons[:, :, 1]

##########################
# Get data to plot state and province boundaries
states_provinces = cfeature.NaturalEarthFeature(category='cultural',
                                                name='admin_1_states_provinces_lakes',
                                                scale='50m',
                                                facecolor='none')

##########################
# Define a 2-D lat/lon index function to find the nearest grid point to a specified coordinate


def lat_lon_2d_index(y, x, lat1, lon1):
    """Calculate the index values of the grid point nearest a given lat/lon point.

    This function calculates the distance from a desired lat/lon point
    to each element of a 2D array of lat/lon values, typically from model output,
    and determines the index value corresponding to the nearest lat/lon grid point.

    y = latitude array
    x = longitude array
    lon1 = longitude point (signle value)
    lat1 = latitude point (single value)

    Returns the index value for nearest lat/lon point on grid

    Equations for variable distiance between longitudes from
    http://andrew.hedges.name/experiments/haversine/, code by Kevin Goebbert.
    """
    r = 6373. * 1000.  # Earth's Radius in meters
    rad = np.pi / 180.
    x1 = np.ones(x.shape) * lon1
    y1 = np.ones(y.shape) * lat1
    dlon = np.abs(x - x1)
    dlat = np.abs(y - y1)
    a = (np.sin(rad * dlat / 2.))**2 + (np.cos(rad * y1) * np.cos(rad * y) *
                                        (np.sin(rad * dlon / 2.))**2)
    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
    d = r * c
    return np.unravel_index(d.argmin(), d.shape)


##############################
# **Insentropic Cross Section**
#
# With the projection set up, we will define the endpoints, extract the
# cross-sectional data, and plot with an inset map.


def cross_section(isentlev, num, left_lat, left_lon, right_lat, right_lon):
    """Plot an isentropic cross-section."""
    # get the coordinates of the endpoints for the cross-section
    left_coord = np.array((float(left_lat), float(left_lon)))
    right_coord = np.array((float(right_lat), float(right_lon)))

    # Calculate data for the inset isentropic map
    isent_anal = mcalc.isentropic_interpolation(float(isentlev) * units.kelvin, lev, tmp,
                                                spech, tmpk_out=True)
    isentprs = isent_anal[0]
    isenttmp = isent_anal[1]
    isentspech = isent_anal[2]
    isentrh = mcalc.relative_humidity_from_specific_humidity(isentspech, isenttmp, isentprs)

    # Find index values for the cross section slice
    iright = lat_lon_2d_index(lat, lon, right_coord[0], right_coord[1])
    ileft = lat_lon_2d_index(lat, lon, left_coord[0], left_coord[1])

    # Get the cross-section slice data
    cross_data = mcalc.extract_cross_section(ileft, iright, lat, lon, tmp, uwnd, vwnd, spech,
                                             num=num)
    cross_lat = cross_data[0]
    cross_lon = cross_data[1]
    cross_t = cross_data[2]
    cross_u = cross_data[3]
    cross_v = cross_data[4]
    cross_spech = cross_data[5]

    # Calculate theta and RH on the cross-section
    cross_theta = mcalc.potential_temperature(lev[:, np.newaxis], cross_t)
    cross_rh = mcalc.relative_humidity_from_specific_humidity(cross_spech, cross_t,
                                                              lev[:, np.newaxis])

    # Create figure for ploting
    fig = plt.figure(1, figsize=(17., 12.))

    # Plot the cross section
    ax1 = plt.axes()
    ax1.set_yscale('symlog')
    ax1.grid()
    cint = np.arange(250, 450, 5)

    # Determine whether to label x-axis with lat or lon values
    if np.abs(left_lon - right_lon) > np.abs(left_lat - right_lat):
        cs = ax1.contour(cross_lon, lev[::-1], cross_theta[::-1, :], cint, colors='tab:red')
        cf = ax1.contourf(cross_lon, lev[::-1], cross_rh[::-1, :], range(10, 106, 5),
                          cmap=plt.cm.gist_earth_r)
        ax1.barbs(cross_lon[4::4], lev, cross_u[:, 4::4], cross_v[:, 4::4], length=6)
        plt.xlabel('Longitude (Degrees East)')
    else:
        cs = ax1.contour(cross_lat[::-1], lev[::-1], cross_theta[::-1, ::-1], cint,
                         colors='tab:red')
        cf = ax1.contourf(cross_lat[::-1], lev[::-1], cross_rh[::-1, ::-1], range(10, 106, 5),
                          cmap=plt.cm.gist_earth_r)
        ax1.barbs(cross_lat[::-4], lev, cross_u[:, ::-4], cross_v[:, ::-4], length=6)
        plt.xlim(cross_lat[0], cross_lat[-1])
        plt.xlabel('Latitude (Degrees North)')

    # Label the cross section axes
    plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7,
               fmt='%i', rightside_up=True, use_clabeltext=True)
    cb = plt.colorbar(cf, orientation='horizontal', extend=max, aspect=65, shrink=0.75,
                      pad=0.06, extendrect='True')
    cb.set_label('Relative Humidity', size='x-large')
    plt.ylabel('Pressure (hPa)')
    ax1.set_yticklabels(np.arange(1000, 50, -50))
    plt.ylim(lev[0], lev[-1])
    plt.yticks(np.arange(1000, 50, -50))

    # Add a title
    plt.title(('NARR Isentropic Cross-Section: ' + str(left_coord[0]) + ' N, ' +
               str(left_coord[1]) + ' E  to ' + str(right_coord[0]) + ' N, ' +
               str(right_coord[1]) + ' E'), loc='left')
    plt.title('VALID: {:s}'.format(str(vtimes[0])), loc='right')

    # Add Inset Map
    ax2 = fig.add_axes([0.125, 0.643, 0.25, 0.25], projection=crs)

    # Coordinates to limit map area
    bounds = [(-122., -75., 25., 50.)]

    # Limit extent of inset map
    ax2.set_extent(*bounds, crs=ccrs.PlateCarree())
    ax2.coastlines('50m', edgecolor='black', linewidth=0.75)
    ax2.add_feature(states_provinces, edgecolor='black', linewidth=0.5)

    # Plot the surface
    clevisent = np.arange(0, 1000, 25)
    cs = ax2.contour(tlons, tlats, isentprs[0, :, :], clevisent,
                     colors='k', linewidths=1.0, linestyles='solid')
    plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7,
               fmt='%i', rightside_up=True, use_clabeltext=True)

    # Plot RH
    cf = ax2.contourf(tlons, tlats, isentrh[0, :, :], range(10, 106, 5),
                      cmap=plt.cm.gist_earth_r)

    # Convert endpoints of cross-section line
    left = crs.transform_point(left_coord[1], left_coord[0], ccrs.PlateCarree())
    right = crs.transform_point(right_coord[1], right_coord[0], ccrs.PlateCarree())

    # Plot the cross section line
    plt.plot([left[0], right[0]], [left[1], right[1]], color='r')
    plt.show()


# Finally, call the plotting function and show the map interactively
interact_manual(cross_section,
                isentlev=widgets.IntSlider(min=280, max=330, step=1, value=296,
                                           description='Isentropic Level'),
                num=widgets.IntSlider(min=20, max=500, step=10, value=80,
                                      description='Sampling Points'),
                left_lat=widgets.FloatSlider(min=0., max=90., step=.5, value=37.,
                                             description='West Lat (N)'),
                left_lon=widgets.FloatSlider(min=-150, max=-45, step=0.5, value=-105.,
                                             description='West Lon (E)'),
                right_lat=widgets.FloatSlider(min=0., max=90., step=.5, value=35.5,
                                              description='East Lat (N)'),
                right_lon=widgets.FloatSlider(min=-150, max=-45, step=.5, value=-65.,
                                              description='East Lon (E)'))


1			# Copyright (c) 2008-2017 MetPy Developers.
2			# Distributed under the terms of the BSD 3-Clause License.
3			# SPDX-License-Identifier: BSD-3-Clause
4			"""
5			Isentropic Cross-section
6			========================
7
8			Make an isentropic cross section with MetPy
9
10			With the function `metpy.calc.extract_cross_section`, data can be extracted from 3-dimensional
11			input along a specified line and plotted in a 2-dimensional cross-section.
12
13			.. image:: ../../_static/cross_section.gif
14
15			"""
16
17			################################
18			import cartopy.crs as ccrs
19			import cartopy.feature as cfeature
20			import ipywidgets as widgets
21			from ipywidgets import interact_manual
22			import matplotlib.pyplot as plt
23			from netCDF4 import Dataset, num2date
24			import numpy as np
25
26			import metpy.calc as mcalc
27			from metpy.cbook import get_test_data
28			from metpy.units import units
29
30			#######################################
31			# Getting the data
32			#
33			# In this example, NARR reanalysis data for 18 UTC 04 April 1987 from the National Centers
34			# for Environmental Information (https://nomads.ncdc.noaa.gov) will be used.
35
36			data = Dataset(get_test_data('narr_example.nc', False))
37
38			##########################
39			print(list(data.variables))
40
41			#############################
42			# We will reduce the dimensionality of the data as it is pulled in to remove an empty time
43			# dimension. Additionally, units are required for input data, so the proper units will also
44			# be attached.
45
46
47			# Assign data to variable names
48			dtime = data.variables['Geopotential_height'].dimensions[0]
49			dlev = data.variables['Geopotential_height'].dimensions[1]
50			lat = data.variables['lat'][:]
51			lon = data.variables['lon'][:]
52			lev = data.variables[dlev][:] * units(data.variables[dlev].units)
53			times = data.variables[dtime]
54			vtimes = num2date(times[:], times.units)
55
56			temps = data.variables['Temperature']
57			tmp = temps[0, :] * units.kelvin
58			uwnd = data.variables['u_wind'][0, :] * units(data.variables['u_wind'].units)
59			vwnd = data.variables['v_wind'][0, :] * units(data.variables['v_wind'].units)
60			hgt = data.variables['Geopotential_height'][0, :] * units.meter
61			spech = (data.variables['Specific_humidity'][0, :] *
62			units(data.variables['Specific_humidity'].units))
63
64
65			#########################
66			# Create the coordinate system and projection for plotting the inset map
67			crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)
68			tlatlons = crs.transform_points(ccrs.PlateCarree(), lon, lat)
69			tlons = tlatlons[:, :, 0]
70			tlats = tlatlons[:, :, 1]
71
72			##########################
73			# Get data to plot state and province boundaries
74			states_provinces = cfeature.NaturalEarthFeature(category='cultural',
75			name='admin_1_states_provinces_lakes',
76			scale='50m',
77			facecolor='none')
78
79			##########################
80			# Define a 2-D lat/lon index function to find the nearest grid point to a specified coordinate
81
82
83			def lat_lon_2d_index(y, x, lat1, lon1):
84			"""Calculate the index values of the grid point nearest a given lat/lon point.
85
86			This function calculates the distance from a desired lat/lon point
87			to each element of a 2D array of lat/lon values, typically from model output,
88			and determines the index value corresponding to the nearest lat/lon grid point.
89
90			y = latitude array
91			x = longitude array
92			lon1 = longitude point (signle value)
93			lat1 = latitude point (single value)
94
95			Returns the index value for nearest lat/lon point on grid
96
97			Equations for variable distiance between longitudes from
98			http://andrew.hedges.name/experiments/haversine/, code by Kevin Goebbert.
99			"""
100			r = 6373. * 1000. # Earth's Radius in meters
101			rad = np.pi / 180.
102			x1 = np.ones(x.shape) * lon1
103			y1 = np.ones(y.shape) * lat1
104			dlon = np.abs(x - x1)
105			dlat = np.abs(y - y1)
106			a = (np.sin(rad * dlat / 2.))*2 + (np.cos(rad y1) * np.cos(rad * y) *
107			(np.sin(rad * dlon / 2.))**2)
108			c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
109			d = r * c
110			return np.unravel_index(d.argmin(), d.shape)
111
112
113			##############################
114			# Insentropic Cross Section
115			#
116			# With the projection set up, we will define the endpoints, extract the
117			# cross-sectional data, and plot with an inset map.
118
119
120			def cross_section(isentlev, num, left_lat, left_lon, right_lat, right_lon):
121			"""Plot an isentropic cross-section."""
122			# get the coordinates of the endpoints for the cross-section
123			left_coord = np.array((float(left_lat), float(left_lon)))
124			right_coord = np.array((float(right_lat), float(right_lon)))
125
126			# Calculate data for the inset isentropic map
127			isent_anal = mcalc.isentropic_interpolation(float(isentlev) * units.kelvin, lev, tmp,
128			spech, tmpk_out=True)
129			isentprs = isent_anal[0]
130			isenttmp = isent_anal[1]
131			isentspech = isent_anal[2]
132			isentrh = mcalc.relative_humidity_from_specific_humidity(isentspech, isenttmp, isentprs)
133
134			# Find index values for the cross section slice
135			iright = lat_lon_2d_index(lat, lon, right_coord[0], right_coord[1])
136			ileft = lat_lon_2d_index(lat, lon, left_coord[0], left_coord[1])
137
138			# Get the cross-section slice data
139			cross_data = mcalc.extract_cross_section(ileft, iright, lat, lon, tmp, uwnd, vwnd, spech,
140			num=num)
141			cross_lat = cross_data[0]
142			cross_lon = cross_data[1]
143			cross_t = cross_data[2]
144			cross_u = cross_data[3]
145			cross_v = cross_data[4]
146			cross_spech = cross_data[5]
147
148			# Calculate theta and RH on the cross-section
149			cross_theta = mcalc.potential_temperature(lev[:, np.newaxis], cross_t)
150			cross_rh = mcalc.relative_humidity_from_specific_humidity(cross_spech, cross_t,
151			lev[:, np.newaxis])
152
153			# Create figure for ploting
154			fig = plt.figure(1, figsize=(17., 12.))
155
156			# Plot the cross section
157			ax1 = plt.axes()
158			ax1.set_yscale('symlog')
159			ax1.grid()
160			cint = np.arange(250, 450, 5)
161
162			# Determine whether to label x-axis with lat or lon values
163			if np.abs(left_lon - right_lon) > np.abs(left_lat - right_lat):
164			cs = ax1.contour(cross_lon, lev[::-1], cross_theta[::-1, :], cint, colors='tab:red')
165			cf = ax1.contourf(cross_lon, lev[::-1], cross_rh[::-1, :], range(10, 106, 5),
166			cmap=plt.cm.gist_earth_r)
167			ax1.barbs(cross_lon[4::4], lev, cross_u[:, 4::4], cross_v[:, 4::4], length=6)
168			plt.xlabel('Longitude (Degrees East)')
169			else:
170			cs = ax1.contour(cross_lat[::-1], lev[::-1], cross_theta[::-1, ::-1], cint,
171			colors='tab:red')
172			cf = ax1.contourf(cross_lat[::-1], lev[::-1], cross_rh[::-1, ::-1], range(10, 106, 5),
173			cmap=plt.cm.gist_earth_r)
174			ax1.barbs(cross_lat[::-4], lev, cross_u[:, ::-4], cross_v[:, ::-4], length=6)
175			plt.xlim(cross_lat[0], cross_lat[-1])
176			plt.xlabel('Latitude (Degrees North)')
177
178			# Label the cross section axes
179			plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7,
180			fmt='%i', rightside_up=True, use_clabeltext=True)
181			cb = plt.colorbar(cf, orientation='horizontal', extend=max, aspect=65, shrink=0.75,
182			pad=0.06, extendrect='True')
183			cb.set_label('Relative Humidity', size='x-large')
184			plt.ylabel('Pressure (hPa)')
185			ax1.set_yticklabels(np.arange(1000, 50, -50))
186			plt.ylim(lev[0], lev[-1])
187			plt.yticks(np.arange(1000, 50, -50))
188
189			# Add a title
190			plt.title(('NARR Isentropic Cross-Section: ' + str(left_coord[0]) + ' N, ' +
191			str(left_coord[1]) + ' E to ' + str(right_coord[0]) + ' N, ' +
192			str(right_coord[1]) + ' E'), loc='left')
193			plt.title('VALID: {:s}'.format(str(vtimes[0])), loc='right')
194
195			# Add Inset Map
196			ax2 = fig.add_axes([0.125, 0.643, 0.25, 0.25], projection=crs)
197
198			# Coordinates to limit map area
199			bounds = [(-122., -75., 25., 50.)]
200
201			# Limit extent of inset map
202			ax2.set_extent(*bounds, crs=ccrs.PlateCarree())
203			ax2.coastlines('50m', edgecolor='black', linewidth=0.75)
204			ax2.add_feature(states_provinces, edgecolor='black', linewidth=0.5)
205
206			# Plot the surface
207			clevisent = np.arange(0, 1000, 25)
208			cs = ax2.contour(tlons, tlats, isentprs[0, :, :], clevisent,
209			colors='k', linewidths=1.0, linestyles='solid')
210			plt.clabel(cs, fontsize=10, inline=1, inline_spacing=7,
211			fmt='%i', rightside_up=True, use_clabeltext=True)
212
213			# Plot RH
214			cf = ax2.contourf(tlons, tlats, isentrh[0, :, :], range(10, 106, 5),
215			cmap=plt.cm.gist_earth_r)
216
217			# Convert endpoints of cross-section line
218			left = crs.transform_point(left_coord[1], left_coord[0], ccrs.PlateCarree())
219			right = crs.transform_point(right_coord[1], right_coord[0], ccrs.PlateCarree())
220
221			# Plot the cross section line
222			plt.plot([left[0], right[0]], [left[1], right[1]], color='r')
223			plt.show()
224
225
226			# Finally, call the plotting function and show the map interactively
227			interact_manual(cross_section,
228			isentlev=widgets.IntSlider(min=280, max=330, step=1, value=296,
229			description='Isentropic Level'),
230			num=widgets.IntSlider(min=20, max=500, step=10, value=80,
231			description='Sampling Points'),
232			left_lat=widgets.FloatSlider(min=0., max=90., step=.5, value=37.,
233			description='West Lat (N)'),
234			left_lon=widgets.FloatSlider(min=-150, max=-45, step=0.5, value=-105.,
235			description='West Lon (E)'),
236			right_lat=widgets.FloatSlider(min=0., max=90., step=.5, value=35.5,
237			description='East Lat (N)'),
238			right_lon=widgets.FloatSlider(min=-150, max=-45, step=.5, value=-65.,
239			description='East Lon (E)'))
240

Unidata / MetPy

Pull Request — master (#518)

cross_section() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like