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.

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

Problems with profile calculations and NUCAPS data #973

Closed
kuciauskas opened this issue Nov 28, 2018 · 5 comments
Closed

Problems with profile calculations and NUCAPS data #973

kuciauskas opened this issue Nov 28, 2018 · 5 comments
Labels
Area: Calc Pertains to calculations Area: Docs Affects documentation Type: Bug Something is not working like it should
Milestone

Comments

@kuciauskas
Copy link

kuciauskas commented Nov 28, 2018

I am using NUCAPS sounding datasets to plot SkewTs. Although I can plot the temperature and dewpt temp profiles, I am having trouble with the advanced sounding aspects as provided within the MetPy examples

  1. calculating CAPE and Convective Inhibition (cin) values
    • I am getting consistently unreasonably high cin values and zeros for CAPE
  2. LCL (temperature) and LCL (pressure) values -
  3. unable to plot the parcel profile as a solid black line
  4. unable to annotating text on the skewT plot using the above values
from netCDF4 import Dataset

# Installed Libraries
from IPython import embed as shell
import numpy as np
import netCDF4 as ncdf
import metpy.calc as mpcalc
from metpy.plots import add_metpy_logo, SkewT
from metpy.units import units
from netCDF4 import Dataset

import matplotlib.pyplot as plt

# work on a selected set of files
for path in glob.glob("NUCAPS*.nc"):
    print(path)
    nc = Dataset(path)

    for index in range(len(nc.dimensions['Number_of_CrIS_FORs'])):
        time = nc.variables['Time'][index]
        tym = time/1000.0
        final_time = datetime.datetime.fromtimestamp(tym).strftime('%Y%m%d.%H%M%S')

#    skew = SkewT(plt.figure(figsize=(6, 7)))
        lat = nc.variables[u'Latitude'][index]
        lon = nc.variables[u'Longitude'][index]
        lon_a = lon

# check if lat and lon are within the range
        if ((lat > 12.0 and lat < 16.0)) and ((lon > -19.0 and lon < -15.0)):
           print 'lat: ', lat, 'and lon: ', lon, ' is within limts'
        else:
           continue

        lon = abs(lon)
        lat = "%.1f" % (float(lat))
        lon = "%.1f" % (float(lon))
        lon_a = "%.1f" % (float(lon_a))

        fig = plt.figure(figsize=(13,8))
        add_metpy_logo(fig, 115, 100)
#        skew = SkewT(plt.figure(figsize=(13,8)))
        skew = SkewT(fig, rotation=45)

        temp_K = nc.variables[u'Temperature'][index] * units.kelvin
        temp_C = temp_K.to('degC')
        mixing = nc.variables[u'H2O_MR'][index]
        press = nc.variables[u'Pressure'][index] * units.hPa

        e = mpcalc.vapor_pressure(press, mixing)
        dewp_C = mpcalc.dewpoint(e)
        dewp_K = dewp_C.to('degK')

# Calculate LCL height and plot as black dot
        lcl_pressure, lcl_temperature = mpcalc.lcl(press[0], temp_C[0], dewp_C[0])
        skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')


# Plot the parcel profile as a black line
        prof = mpcalc.parcel_profile(press, temp_K[0], dewp_K[0]).to('degC')
        prof_K = mpcalc.parcel_profile(press, temp_K[0], dewp_K[0])
        skew.plot(press, prof, 'k', linewidth=4)

# plot out the skew-T paramters
# Add the relevant special lines
#        skew.plot(press, temp.to('degC'), color='tab:red', linewidth=3, alpha=0.5)
        skew.plot(press, temp_C, color='tab:red', linewidth=3, alpha=0.5)
#        skew.plot(press, dewp.to('degC'), color='tab:blue', linewidth=3, alpha=0.5)
        skew.plot(press, dewp_C, color='tab:green', linewidth=3, alpha=0.5)
        skew.ax.set_xlim(-50, 50)
        skew.ax.set_ylim(1000, 100)
        skew.plot_dry_adiabats()
        skew.plot_moist_adiabats()

# An example of a slanted line at constant T -- in this case the 0
# isotherm
        skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)

#        skew.plot_mixing_lines()
# Shade areas of CAPE and CIN
#        skew.shade_cin(press, temp, prof)
        skew.shade_cape(press, temp_K, prof)

# calculate Convective inhibition (CIN) and Convective available potential energy (CAPE)
#        cin = mpcalc.cape_cin(press, temp, dewp, prof)
        cin = mpcalc.cape_cin(press, temp_C, dewp_C, prof_K)
        lfc = mpcalc.lfc(press, temp_K, dewp_K, prof)
        cin_val = cin[1].magnitude
        cin_val = "%.2f" % (float(cin_val))
        print('cin: ',cin_val, 'lfc: ', lfc)

        plt.ylabel('log P')
        title = "DTG: ", final_time, "    Lat: ", lat, "Lon: ", lon_a, "CI: ", cin_val
        plt.title(title, fontsize=12)
        index = "%03d" % (int(index+1))

#        plt.savefig('test_{}.png'.format(index))
        temp_fn = "NUCAPS_"+str(final_time)+".lat"+str(lat)+"lon"+str(lon)
        plt.savefig('{}.png'.format(temp_fn))
        plt.close()
    plt.show()
@dopplershift
Copy link
Member

Sample Image:
nucaps_20181118 192304 lat14 5lon18 8

And data.

@dopplershift dopplershift changed the title SkewT: dataset, plot, and code Problems with profile calculations and NUCAPS data Nov 28, 2018
@dopplershift dopplershift added Type: Bug Something is not working like it should Area: Calc Pertains to calculations labels Nov 28, 2018
@dopplershift
Copy link
Member

So, I've found two issues here.

First, the NUCAPS data comes out in pressure-order, meaning low pressure to high pressure. MetPy's profile calculation functions, like CAPE and parcel_profile, assume data start at the bottom of the atmosphere and go upwards. For now, you can work around this by flipping the profiles (using slicing with [::-1]) when you read them in:

        temp_K = nc.variables[u'Temperature'][index][::-1] * units.kelvin
        mixing = nc.variables[u'H2O_MR'][index][::-1]
        press = nc.variables[u'Pressure'][index][::-1] * units.hPa

Ideally, we could handle some of this automatically, but it starts making calculations convoluted. We could definitely document this more clearly. We may also consider a helper function that checks the data once and adjusts as necessary.

Second, netCDF4-python is returning masked arrays for the data (even though there are no masked points). For some reason, that needs investigating, the masked arrays are breaking the CAPE/CIN calculation. As a workaround, you can adjust your code to open the file as follows:

    nc = Dataset(path)
    nc.set_always_mask(False)

That has things working sensibly for me now:
sounding

@dopplershift dopplershift added the Area: Docs Affects documentation label Nov 28, 2018
@dopplershift dopplershift added this to the 0.10 milestone Nov 28, 2018
@kuciauskas
Copy link
Author

kuciauskas commented Nov 28, 2018 via email

@kuciauskas
Copy link
Author

Ryan:

In your comments above:
instead of: nc.set_always_mask(False) <I don't see that attribute for nc>
shouldn't it be: nc.set_auto_mask(False)

Additionally, the cin values after executing:
cin = mpcalc.cape_cin(press, temp_C, dewp_C, prof)
are NaN for the 1st index#, on the order of -26 for the 2nd index

@dopplershift
Copy link
Member

What version of netcdf4 do you have? (Try: conda list netcdf4) The latest is 1.4.2. Then again, if you're running with an older netcdf4, you might not have the problem with it always returning a masked array.

Regarding the cin values, you're not quite using it right. That should be:

cape, cin = mpcalc.cape_cin(press, temp_C, dewp_C, prof)

the function calculates both CAPE and CIN and returns them both. So the result you're seeing is a NaN for CAPE, because there is no positive area, and -26 for CIN because that's the area below (I think) the LCL.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Calc Pertains to calculations Area: Docs Affects documentation Type: Bug Something is not working like it should
Projects
None yet
Development

No branches or pull requests

3 participants