The excerpt from the text forgets to mention that you assume Local Thermodynamic Equilibrium, and not full Thermodynamic Equilibrium, so to make it possible to define point to point (or from region to region) an EoS.
If there is no sense of being 'close' to thermodynamical equilibrium, it is simply impossible to talk about EoS, pressure and the like from the "Hydrodynamics as Local Thermal Equilibrium".
From the strict thermodynamic view, you can't talk of anything time-dependent. All piston an thermal cycles use the abuse of talking of 'quasi-equilibrium' without really defining it, and simply postulate that you can use full traditional thermodynamics all around. From a 'rigorous' point of view, the only thing that you can talk in thermodynamics are stationary, homogeneous systems, which are in the 'thermodynamical limit' ( $N,V,S...\rightarrow \infty$ but $N/V, S/V...$ fixed), so no time dependence.
The idea of flow is that even if the fluid flows very fast in relation to a fixed observer, if you go to the rest frame of that piece of fluid, you can talk of thermodynamic equilibrium close to that small part of the fluid.
I believe that the best way to understand how this all works is through Boltzmann Equation, which I can develop latter if you wish so.
Edit(Complementing the answer, as asked):
So, you have boltzmann equation:
$
\frac{\partial f}{\partial t}+\vec v \cdot \frac{\partial f}{\partial \vec r}+\vec F\cdot \frac{\partial f}{\partial \vec p} = \int d^3 \vec p_0 d\Omega\ g\ \sigma(g,\Omega) (f'f_1' - ff_1)
$
Where:$g=|\vec p - \vec p_0|$, $\sigma(g,\Omega)$ is the differential cross section between gas molecules, and the primed distributions are evaluated with the momentum that corresponds to an out-going (solid)-angle $\Omega$, with ingoing momenta $\vec p$ and $\vec p_0$. The normalization is $\int d^3\vec r\ d^3 \vec p\ f(t,\vec r,\vec p)=N$, where N is the total number of particles.
We believe that this equation provides a good description to 1-particle distribution function, in phase space, of a a rarefied gas composed with hard-spheres(i.e. hard, short range, repulsive potential, with only elastic collisions). Putting aside whether it's justified or not to model a gas this way, simply believe for the moment that it works.
Now you want to model a gas inside a box as beeing in full thermodynamical equilibrium. Equilibrium is when you have stationary, homogeneous material. So, you want to look for solutions of Boltzmann equation that have this kind of symmetry, and thus:
$f(t,\vec r, \vec p) = \frac{N}{V}Id_V(\vec r)f_0(\vec p)$
So, the most of inhomogeneity that there may be is an indicator function that says that outside the box, there is no gas. We are also supposing that the only external forces are on the walls of the box, and so, in the bulk of the gas we have $\vec F=0$
Now we feed this ansatz to the Boltzmann equation and see what happens. Now, from the above assumptions $\partial f/\partial t=0$, $\partial f/\partial \vec r = 0$ and $\vec F=0$ on the bulk. This gives us:
$
\frac{\partial f}{\partial t}+\vec v \cdot \frac{\partial f}{\partial \vec r}+\vec F\cdot \frac{\partial f}{\partial \vec p} = 0 = \int d^3 \vec p_0 d\Omega\ g\ \sigma(g,\Omega) (f'f_1' - ff_1)
$
So we need to kill the collision kernel in order to satisfy the Boltzmann equation. The easiest way is to nullify the subtraction inside it by putting:
$
f_0(\vec p)f_0(\vec p_1)=f_0(\vec p')f_0(\vec p_1')
$
For all possible (elastic) collision outcomes. Now comes the the smart point. Lets take the $\log$ of the above expression.
$
\log f_0(\vec p) + \log f_0(\vec p_1)=\log f_0(\vec p') + \log f_0(\vec p_1')
$
If $\log f_0$ is function only of additive conserved quantities on the collision, we get the relation above for free! (Ok, not completely for free, its possible to show that this is essentially the only way to do it)
Now, for elastic binary collisions, we have only 3 conserved quantitites: Mass, Linear Momentum (because we believe that there is no relevant rotation) and kinectic energy.
Now we write:
$\log f_0(\vec p) = A\frac{\vec p^2}{2m}+ \vec B \cdot \vec p + Cm$
Massaging the above expression and we use integrability conditions, we may write:
$\log f_0(\vec p) = -\frac{ (\vec p-\vec p_0)^2}{2m\sigma^2}+ \log N_0$
In the case of the box, we know that the box isn't moving (equivalently, it's locally isotropic), so we put $\vec p_0=0$ and we get the Boltzmann Distribution as a solution to Boltzmann equation for equilibrium conditions. Further, we can identify $\sigma^2 = k_BT$ and we close the identification.
Now, to hydrodynamics. To find hydrodynamical equation from Botlzmann equation, we "Take Moments" from it, i.e., e multiply it by powers of the linear momentum, and integrate in momentum, so we get equations for things that live in usual 3D space.
Multiplying by $\chi(\vec p)$ and integrating:
$
\frac{\partial}{\partial t}\left(\int d^3\vec p\ \chi(\vec p) f\right) + \frac{1}{m} \nabla_{\vec r} \cdot \left(\int d^3\vec p\ \chi(\vec p)\vec p f\right) + \vec F \cdot \left(\int d^3\vec p\ \chi(\vec p) \frac{\partial f}{\partial \vec p}\right) = \int d^3 \vec p\ d^3 \vec p_0\ d\Omega\ g\chi(\vec p)\ \sigma(g,\Omega) (f'f_1' - ff_1)
$
It's possible to show that if $\chi(\vec p)$ is a conserved quantity on binary collisions, the last term is $=0$, so that's what we are going to look for. Choosing $\chi(\vec p)=m$, we arrive at:
$
\frac{\partial}{\partial t}\left(\int d^3\vec p\ m f\right) + \nabla_{\vec r} \cdot \left(\int d^3\vec p\ \vec p f\right) = 0
$
Identifying $\rho = \int d^3\vec p\ m f$ as the mass density and $\int d^3\vec p\ \vec p f = \vec j = \rho <\vec v> = \rho \vec u$ the mass current, we have the continuity equation for mass density:
$
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec v) = 0
$
Setting $\chi(\vec p) = \vec p$ we arrive at:
$
\frac{\partial }{\partial t}(\rho \vec u) + \nabla \cdot \Sigma - \rho F = 0
$
Where $\Sigma_{ij} = \int d^3\vec p\ p_i p_j f(t,\vec r,\vec p)$. Now we can decompose $\vec p = m<\vec v> + \delta \vec p = m\vec u + \delta \vec p$, where we identify the average velocity as $\vec u = \frac{1}{\rho} \vec j$. This average velocity is what we identify as the fluid velocity. Going back to the last equation we have:
$
\frac{\partial }{\partial t}(\rho \vec u) + \nabla \cdot \left(\rho \vec u \otimes \vec u + \Pi
\right) = \vec f
$
Where, finally, we identify $\Pi_{ij} = \int d^3\vec p\ \delta p_i \delta p_j f(t,\vec r,\vec p)$ as the stress tensor. The convective par is already there, and we now can appreciate the conection between kinectic theory and hydrodynamics. Comming back to Boltzmann Distribuition:
$f=n(t,\vec r)f_0(\vec p)$
$f_0(\vec p) = \frac{1}{(2\pi m k_BT)^{3/2}}e^{-\frac{ (\vec p-\vec p_0)^2}{2mk_BT}}$
We said that for Thermodynamics, we had $n$,$T$ and $\vec p_0$ constant all along the gas. For Hydrodynamics, we try to retain that functional form, and relax this assumption, i.e., we try (again) to find solutions to Boltzmann equation with the above form, but with $T(t,\vec r)$ and $\vec p(t,\vec r)$ possibly have some dependence in time and space, and so we talk about local thermal equilibrium, since we try to keep, locally, an equilibrium distribution.
If we do that, we end up with $\rho = m\times n(t,\vec r)$ and $\vec u = \frac{\vec p_0}{m}$, which wasn't totally unexpected, and $\vec p - m\vec u= \delta \vec p$, so the Boltzmann distribution measures the (local) fluctuation of velocity. Now computing the stress tensor:
$\Pi_{ij} = \frac{n}{(2\pi m k_BT)^{3/2}}\int d^3\vec p\ \delta p_i \delta p_j e^{-\frac{\delta \vec p^2}{2mk_BT}}$
It' not too difficult to see that the stress tensor above is proportional to the identity tensor, and we identify $\Pi_{ij} = p \delta_{ij}$, and since we have a relation between pressure, density and temperature, we have an EoS. If you plug that on the original equation with $\Pi$, you end up with Euler equation for Fluid Dynamics. So you can think that Euler equation is an evolution equation for something that is in strict local equilibrium.
Also, if you look closely, you will see that the probability distribution only care about the velocity fluctuation $\delta \vec p/m$, and not the actual velocity of the fluid $\vec u$. Here enters your question about the fluid flow:
But in fluid dynamics fluids are flowing and that too rapidly and from intuition there will not be any thermodynamic equilibrium during fluid flow.
From the fluid standpoint, the average velocity is not important to the thermodynamics, only the fluctuations around this average velocity.
Chemical equilibrium is not being considered here, since we are supposing that the fluid have a single chemical species, so it's naturally in chemical equilibrium.
Now, beyond Euler equation:
One very (strong) assumption that we made was that the fluid had the distribution in phase space that was locally Maxwell-Boltzmann. What would happen if we dropped this assumption?
Generally, we can't solve (or can only solve numerically) Boltzmann Equation except on very special cases, so, as any good physicist, we go to the next best thing: Approximate Solutions
What happens if our system is not on equilibrium but close to equilibrium? It should be possible to write $f=f_0\phi$ where $\phi \approx 1$. Now, you would like to find some parameter that you could use to to some kind of perturbation expansion around it. This parameter is essentially the Knudsen Number of the system. If you do this, essentially, the only thing that should change here is the stress tensor, which depends explicitly of the form of the distribution on phase space.
The Knudsen Number is essencially a measure of how far the "microscopic" scale of your system is far from the "macroscopic" scale. If they are sufficiently far apart, i.e. $Kn << 1$, an macroscopic, or Hydrodynamical, description of your system should be good.
The zero-th order on Kn should be $f_0$, so you seek something like $\phi = 1+Kn \phi_1 + (Kn)^2 \phi_2 + ...$
You can carry out this calculation, which is rather lengthy, and what you find (if I remember correctly) is that in first order on the Knudsen Number, you find the Navier-Stokes Stress tensor, which in turn bring you to Navier-Stokes equation, with the bulk and shear viscosity coefficients.
Not only this, you can calculate the dependence on the density and temperature for this coeficients, so you not only have the general form of the evolution equation, but you also have an "EoS" in the extended sense, so to encompass also the viscosity coeffients.
So, the idea of using this method is to define pressure, temperature and the like on the "Equilibrium" part of the distribution, and viscosity and any other kind of effect on the "Non-Equilibrium" part. In this sense you can talk about thermodynamics since you are near (local) equilibrium, even you are not exactly on equilibrium.
So, this is one way to see hydrodynamics as 'mean kinetic theory', and also as an (almost) local thermodynamics. There are also other ways to do it. One is to study Non-Equilibrium Thermodynamics as a macroscopic (in the same sense as classical thermodynamics) mean theory. This is done, in the linear theory, by De Groot & Mazur.
I hope that I have clarified some of your questions. I believe this is a very interesting subject, and I like it very much myself.
