Skip to content

Commit

Permalink
fix overlay figure
Browse files Browse the repository at this point in the history
  • Loading branch information
observingClouds committed Oct 9, 2023
1 parent 57ef62d commit 1f3671f
Showing 1 changed file with 26 additions and 21 deletions.
47 changes: 26 additions & 21 deletions scripts/plot_radiosonde_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -74,39 +74,44 @@ def get_args():

def data4skewt(direction, data):
if direction == "ascent":
p = data["p"].isel({"sounding": 0}).pint.to(units.Pa).values * units.Pa
T = data["ta"].isel({"sounding": 0}).pint.to(units.K).values * units.K
Td = data["dp"].isel({"sounding": 0}).pint.to(units.K).values * units.K
p = data["p"].isel({"sounding": 0}).metpy.to(units.Pa).values * units.Pa
T = data["ta"].isel({"sounding": 0}).metpy.to(units.K).values * units.K
Td = data["dp"].isel({"sounding": 0}).metpy.to(units.K).values * units.K
wind_speed = (
data["wspd"]
.isel({"sounding": 0})
.pint.to(units.meter / units.second)
.metpy.to(units.meter / units.second)
.values
* units.meter
/ units.second
)
wind_dir = (
data["wdir"].isel({"sounding": 0}).pint.to(units.degrees).values
data["wdir"].isel({"sounding": 0}).metpy.to(units.degrees).values
* units.degrees
)
u, v = mpcalc.wind_components(wind_speed, wind_dir)

elif direction == "descent":
p = np.flip(data["p"].isel({"sounding": 0}).pint.to(units.Pa).values) * units.Pa
T = np.flip(data["ta"].isel({"sounding": 0}).pint.to(units.K).values) * units.K
Td = np.flip(data["dp"].isel({"sounding": 0}).pint.to(units.K).values) * units.K
p = (
np.flip(data["p"].isel({"sounding": 0}).metpy.to(units.Pa).values)
* units.Pa
)
T = np.flip(data["ta"].isel({"sounding": 0}).metpy.to(units.K).values) * units.K
Td = (
np.flip(data["dp"].isel({"sounding": 0}).metpy.to(units.K).values) * units.K
)
wind_speed = (
np.flip(
data["wspd"]
.isel({"sounding": 0})
.pint.to(units.meter / units.second)
.metpy.to(units.meter / units.second)
.values
)
* units.meter
/ units.second
)
wind_dir = (
np.flip(data["wdir"].isel({"sounding": 0}).pint.to(units.degrees).values)
np.flip(data["wdir"].isel({"sounding": 0}).metpy.to(units.degrees).values)
* units.degrees
)
u, v = mpcalc.wind_components(wind_speed, wind_dir)
Expand Down Expand Up @@ -148,11 +153,11 @@ def main():
launch_time = np.datetime64(data.launch_time.values[0], "m")
timestamp = str(np.datetime64(launch_time, "m")).replace(":", "")

lats = data.isel(sounding=snd).lat.pint.to("degree_north")
lons = data.isel(sounding=snd).lon.pint.to("degree_east")
lats = data.isel(sounding=snd).lat.metpy.to("degree_north")
lons = data.isel(sounding=snd).lon.metpy.to("degree_east")
alts = (
(data.isel(sounding=snd)[args["altitude"]])
.pint.to(units.km)
.metpy.to(units.km)
.pint.dequantify()
.values
)
Expand All @@ -172,7 +177,7 @@ def main():
alpha=0.5,
linestyle="--",
)
traj = plt.scatter(
traj = ax.scatter(
lons,
lats,
c=alts,
Expand Down Expand Up @@ -211,32 +216,32 @@ def main():
)

ax[0].plot(
data.isel(sounding=snd).ta.pint.to(units.K),
data.isel(sounding=snd).ta.metpy.to(units.K),
alts,
color=quantities_color,
)
ax[1].plot(
data.isel(sounding=snd).dp.pint.to(units.K),
data.isel(sounding=snd).dp.metpy.to(units.K),
alts,
color=quantities_color,
)
ax[2].plot(
data.isel(sounding=snd).rh.pint.to(units.percent),
data.isel(sounding=snd).rh.metpy.to(units.percent),
alts,
color=quantities_color,
)
ax[3].plot(
data.isel(sounding=snd).mr.pint.to(units.percent),
data.isel(sounding=snd).mr.metpy.to(units.percent),
alts,
color=quantities_color,
)
ax[4].plot(
data.isel(sounding=snd).wspd.pint.to(units.meter / units.second),
data.isel(sounding=snd).wspd.metpy.to(units.meter / units.second),
alts,
color=quantities_color,
)
ax[5].plot(
data.isel(sounding=snd).wdir.pint.to(units.degrees),
data.isel(sounding=snd).wdir.metpy.to(units.degrees),
alts,
color=quantities_color,
)
Expand Down Expand Up @@ -306,7 +311,7 @@ def find_nearest(array, value):
# Search for levels by providing pressures
# (levels is the coordinate not pressure)
pres_vals = (
data["p"].isel({"sounding": 0}).pint.to(units.hPa)
data["p"].isel({"sounding": 0}).metpy.to(units.hPa)
).values # Convertion from Pa to hPa

closest_pressure_levels = np.unique(
Expand Down

0 comments on commit 1f3671f

Please sign in to comment.