2

I am trying to use a metpy.calc function however every time I attempt to assign units using metpy.units, it will not run. It will sit there and continue running for hours if I let it. The metpy function has worked perfectly fine for me in the past with the same dataset, but for some reason will no longer work. I have tried updating metpy and spyder and have also used multiple laptops to ensure it wasn't just occurring on one. I also tried the metpy.units two different ways. It runs using units('kg/kg') however won't actually assign the units to the array. Any thoughts on how to get this to work?

This is my code:

import numpy as np 
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import metpy
import metpy.calc as mpcalc
from metpy.units import units

f = Dataset('C:/Users/hkenward/Desktop/MS/cm1out.nc', mode = 'r', format="NETCDF4")
mrv=f.variables['qv'][:]
pt=f.variables['th'][:]
p=f.variables['prs'][:]
p=p/100.
mrv=mrv*units('kg/kg')
p=p * units.hPa
pt=pt*units.kelvin

e=mpcalc.vapor_pressure(p,mrv)
Td=mpcalc.dewpoint(e)
T=mpcalc.temperature_from_potential_temperature(p, pt)
Te=mpcalc.equivalent_potential_temperature(p, T, Td)
DopplerShift
  • 5,472
  • 1
  • 21
  • 20
Holley Kenward
  • 21
  • 2
  • What does "completely breaks my script" mean here? What's your code? What's the error you're getting? – AKX May 19 '20 at 15:51
  • Can you post the code that's failing for you? That's going to be necessary to figure out what's going wrong. https://stackoverflow.com/help/minimal-reproducible-example – DopplerShift May 19 '20 at 18:23

2 Answers2

2

This is a compatibility issue with the Pint units framework (used by MetPy) and masked arrays (which is what netCDF4-python spits out by default). The easiest work around is to multiply your units on the left like:

pt = units.kelvin * pt

That makes everything work. Hopefully things like this will get better in the future by using things like Xarray.

DopplerShift
  • 5,472
  • 1
  • 21
  • 20
1

metpy.calc and metpy.units work fine for me (here on Repl.it):

import metpy.calc as mc
from metpy.units import units

print(mc.density(
  pressure=(1500 * units.pascal),
  temperature=(26 * units.celsius),
  mixing=0.5,
))

outputs

0.014524968314683053 kilogram / meter ** 3
AKX
  • 152,115
  • 15
  • 115
  • 172