Your trajectories must contain rewards, so I'm assuming you've forgotten them in your original post, i.e., we must have $$\tau_j = (s_0^j, a_0^j, r_1^j, ..., s_{N_j}, a_{N_j}, r_{N_j+1})$$
Given that you have access to the full trajectories, I would use the Monte Carlo estimates. You want to use TD methods when you need to estimate $Q^\pi$ incrementally as new transitions $(s_t, a_t, r_{t+1})$ arrive.
You can estimate $Q^\pi$ as follows (adapted from the Sutton&Barto book, chapter 5.1):
Initialise $\text{Returns}(s,a)$ to an empty list for all $(s, a)$ pairs in the trajectories.
For each $j$:
Initialise the cumulative discounted reward $G$ to 0.
For $t$ in ${N_j, N_{j-1}, ..., 0}$:
Set $G = \gamma G + r_{t+1}$
If $(s_t, a_t)$ is not in $((s_0, s_0), ... (s_{t-1}, a_{t-1}))$:
Append $G$ to $\text{Returns}(s_t, a_t)$.
Set $Q^\pi(s_t, a_t) = \text{Average}(\text{Returns}(s_t, a_t))$
For completeness the method is called "First visit Monte Carlo" because we only update the estimate of $Q^\pi$ using first visits to an $(s, a)$ pair.