There are a few issues to unpack here.
First, note that what you have defined is not really the charge density. Let's write out this quantity explicitly:
\begin{equation}
\hat{\rho}(x',y',z', \theta) = \lambda \delta(x'-x(\theta)) \delta(y'-y(\theta)) \delta(z'-z(\theta)),
\end{equation}
Note that I have defined $\hat\rho$ with an explicit dependence on $\theta$. This isn't actually really the charge density, $\rho$. A fixed value of $\theta$ is a single mathematical point on the curve. If the total charge of the whole curve is finite, then the charge at any one infinitesimally small point is actually zero. The total charge density should really "sum" the charge densities from every point on the curve. I put "sum" in scare quotes because the curve is continuous, so we really should integrate. The actual charge density (assuming the charge is uniformly distributed with $\theta$ -- see the footnote below for more on this) is$^\star$
\begin{equation}
\rho(x', y', z') = \int_0^{2\pi} {\rm d} \theta \hat{\rho}(x', y', z', \theta) = \lambda \int_0^{2\pi} {\rm d} \theta \delta(x'-x(\theta)) \delta(y'-y(\theta)) \delta(z'-z(\theta))
\end{equation}
Now, let's apply a little dimensional analysis to our expression for the charge density. $\rho$ needs to have dimensions of charge density (charge per unit volume). Meanwhile, $\delta(\alpha)$ has the same dimensions as $1/\alpha$. Since $x$, $y$, and $z$ all have units of length, the product $\delta(x'-x(\theta)) \delta(y'-y(\theta)) \delta(z'-z(\theta))$ has units of $1/({\rm length})^3$, or $1/{\rm volume}$. Note I've implicitly defined the functions $x(\theta)$, $y(\theta)$, and $z(\theta)$ to be the formulas given in your question; e.g. $z(\theta)= -\sin 3\theta$. Finally, ${\rm d} \theta$ has the same dimensions as $\theta$. Therefore, for dimensional consistency, $\lambda$ must have dimensions of charge per unit $\theta$. In general, you will need to worry about how you choose to parameterize your curve in order to fix $\lambda$. In this case, clearly $0\leq \theta < 2\pi$ is dimensionless. Therefore, $\lambda$ has units of charge.
Now we can think about how to fix $\lambda$. The total charge of a system, which I will denote as $Q_{\rm tot}$, is the volume integral of the charge density
\begin{equation}
Q_{\rm tot} \equiv \int {\rm d} x' {\rm d} y' {\rm d} z' \rho(x', y', z')
\end{equation}
Incidentally, it does not matter what symbol we use to denote this quantity. We can use $Q_{\rm tot}$, $Q$, or $q$, or any other symbol, so long as we are clear on the meaning. There is no intrinsic difference between the notation $Q$ and $q$. I will stick with $Q_{\rm tot}$ in this answer, however, for clarity.
Putting this all together,
\begin{eqnarray}
Q_{\rm tot} &=& \lambda \int_0^{2\pi} {\rm d} \theta \int {\rm d} x' {\rm d} y' {\rm d} z' \delta(x'-x(\theta)) \delta(y'-y(\theta)) \delta (z'-z(\theta)) \\
&=& \lambda \int_0^{2\pi} {\rm d} \theta\\
&=& 2\pi \lambda
\end{eqnarray}
Therefore
\begin{equation}
\lambda = \frac{Q_{\rm tot}}{2\pi}
\end{equation}
$^\star$ This step hides a very important subtlety, since it's not completely obvious what measure factor to use in the integral over $\theta$. Ultimately this requires some physical input -- how is the charge distributed on the trefoil? I'm assuming here that the charge density is uniform as a function of $\theta$, based on the way you phrased the question. However, if the charge density is uniform as a function of length, then you need to include a Jacobian factor converting from arc length to $\theta$ in this step, and carry it through consistently. Depending on how you define things, the dimensions of $\lambda$ may be charge per length, rather than charge, in this case.
In a bit more detail, first define the arc length of the curve as a function of $\theta$
\begin{equation}
L(\theta) = \int_0^\theta {\rm d}\vartheta \sqrt{\dot{x}(\vartheta)^2 + \dot{y}(\vartheta)^2 + \dot{z}(\vartheta)^2}
\end{equation}
where $\dot{x}(\theta)\equiv {\rm d}x/{\rm d}\theta$. The total length of the curve is $L_t = L(2\pi)$.
Then instead of the expression for $\rho$ above, we want to use
\begin{equation}
\rho(x',y',z') = \lambda \int_0^{2\pi} {\rm d} \theta \left|\frac{{\rm d}L}{{\rm d} \theta}\right| \hat{\rho}(x',y',z',\theta)
\end{equation}
Equivalently we could define $\hat{\rho}$ so it was a function of $L$ instead of $\theta$ (e.g. parameterize the curve by its arc length), and then define $\rho$ as an integral over $L$ (with no Jacobian factor) instead of an integral over $\theta$ (with a Jacobian factor)
\begin{equation}
\rho(x',y',z') = \lambda \int_0^{L_t} {\rm d} L \hat{\rho}(x',y',z',L)
\end{equation}
After this, one can follow the logic given in the body of the answer. I believe after turning the crank, the net result should be that $\lambda=Q_{\rm tot}/L_t$.
I believe (though I didn't check explicitly) if you pick, say, 100 uniformly spaced values of $\theta$ and plot the $x,y,z$ coordinates, you'll find that the points tend to "bunch up" near the "curvy parts" of the trefoil, and "spread out" on the "straight parts." On the other hand, if you pick 100 uniformly spaced values of $L(\theta)$ and plot the $x,y,z$ coordinates, the points will be uniformly distributed along the curve with no bunching or spreading out (by construction). Intuitively, the reason you need to worry about this Jacobian factor is to account for the difference in these discretizations.
