def METHOD_NAME(self, weatherFile, defNul=None):
    """
    Read weather variables from a netCDF file
    """
    if defNul is None:
        defNul = np.nan
    # TODO: it'd be cool to use some kind of units package
    # TODO: extract partial pressure directly (q?)
    with netcdf.netcdf_file(weatherFile, 'r', maskandscale=True) as f:
        spvar = f.variables['P_PL']
        temp = f.variables['T_PL']
        humid = f.variables['RH_PL']
        geohvar = f.variables['GHT_PL']
        lon0 = f.STAND_LON.copy()
        lat0 = f.MOAD_CEN_LAT.copy()
        lat1 = f.TRUELAT1.copy()
        lat2 = f.TRUELAT2.copy()
        checkUnits(spvar.units.decode('utf-8'), 'pressure')
        checkUnits(temp.units.decode('utf-8'), 'temperature')
        checkUnits(humid.units.decode('utf-8'), 'relative humidity')
        checkUnits(geohvar.units.decode('utf-8'), 'geopotential')
        # _FillValue is not always set, but when it is we want to read it
        tNull = getNullValue(temp)
        hNull = getNullValue(humid)
        gNull = getNullValue(geohvar)
        pNull = getNullValue(spvar)
        sp = spvar[0].copy()
        temps = temp[0].copy()
        humids = humid[0].copy()
        geoh = geohvar[0].copy()
        spvar = None
        temp = None
        humid = None
        geohvar = None
    # Projection
    # See http://www.pkrc.net/wrf-lambert.html
    earthRadius = 6370e3  # <- note Ray had a bug here
    p1 = CRS(proj='lcc', lat_1=lat1,
             lat_2=lat2, lat_0=lat0,
             lon_0=lon0, a=earthRadius, b=earthRadius,
             towgs84=(0, 0, 0), no_defs=True)
    self._proj = p1
    temps[temps == tNull] = np.nan
    sp[sp == pNull] = np.nan
    humids[humids == hNull] = np.nan
    geoh[geoh == gNull] = np.nan
    self._t = temps
    self._rh = humids
    # Zs are problematic because any z below the topography is nan.
    # For a temporary fix, I will assign any nan value to equal the
    # nanmean of that level.
    zmeans = np.nanmean(geoh, axis=(1, 2))
    nz, ny, nx = geoh.shape
    Zmeans = np.tile(zmeans, (nx, ny, 1))
    Zmeans = Zmeans.T
    ix = np.isnan(geoh)
    geoh[ix] = Zmeans[ix]
    self._zs = geoh
    if len(sp.shape) == 1:
        self._p = np.broadcast_to(
            sp[:, np.newaxis, np.newaxis], self._zs.shape)
    else:
        self._p = sp