In REMD, multiple replicas of the system are simulated "simultaneously" at various temperatures. Using a (Metropolis) Monte Carlo scheme, the replicas 'exchange' temperatures; this 'exchange' scheme will effectively ensure that  the configurations sampled at a particular temperature belong to the equilibrium ensemble at that corresponding temperature. So, your implication in the question, that configuration sampled each temperature  are NOT at equlibrium/boltzmann-weighted is _incorrect_.  (See Sec-14.1 page 319, first paragraph of the book "Understanding Molecular Simulation: From Algorithms to Applications" by Frenkel and Smit, 2nd Edition).

Prabhakar's response is correct in that the sub-ensembles you can extract by looking at all snapshots from a single temperature (or state if you are not doing T-REMD) are valid canonical ensembles by virtue of having used Metropolized Monte Carlo to evaluate replica exchanges.

However, you can actually do a little better than this by reweighting snapshots from all replicas to give improved statistics at each temperature (or state).  As the replica temperatures increase, their weight in the lowest temperature replica will decrease rapidly, however, so the benefit is limited.

You can use either WHAM or MBAR to do this reweighting.  If I was doing it, I would use MBAR (you can find a Python library for doing MBAR calculations on Github: https://github.com/choderalab/pymbar and some common examples there as well for dealing with parallel tempering simulations (T-REMD): https://github.com/choderalab/pymbar/tree/master/examples/parallel-tempering-2dpmf

Ergodicity theorem states that if sampling sufficiently long (how long it is very difficult to say) you will end up with averaging in both in ensemble and in time that would correspond to experiment averages, but if trajectories are not very long it might need some averaging maybe not Boltzmann but a partial Student additional averaging at each temperature. However it is always better to run longer and more accurate runs, since additional averaging can mess up even further everything and introduce artifacts. Try to chose time constants for pressure and temperature that their ratio would give irrational number and pressure coupling should be at least 3 times larger (for each replica and then try to average all of them instead of using mathematical tweaking). Additionally you have to not use Berendsen method (which is more stable with such time constants, but is lacking proper sampling).