Multispecies kinematic flow models are defined by systems of N strongly coupled, nonlinear first-order conservation laws, where the solution is a vector of N partial volume fractions or densities. These models arise in various applications including multiclass vehicular traffic and sedimentation of polydisperse suspensions. The solution vector should take values in a set of physically relevant values (i.e., the components are nonnegative and sum up at most to a given maximum value). It is demonstrated that this set, the so-called invariant region, is preserved by numerical solutions produced by a new family of high-order finite volume numerical schemes adapted to this class of models. To achieve this property, and motivated by [X. Zhang, C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, J. Comput. Phys. 229 (2010) 3091-3120], a pair of linear scaling limiters is applied to a high-order weighted essentially non-oscillatory (WENO) polynomial reconstruction to obtain invariant-region-preserving (IRP) high-order polynomial reconstructions. These reconstructions are combined with a local Lax-Friedrichs (LLF) or Harten-Lax-van Leer (HLL) numerical flux to obtain a high-order numerical scheme for the system of conservation laws. It is proved that this scheme satisfies an IRP property under a suitable Courant-Friedrichs-Lewy (CFL) condition. The theoretical analysis is corroborated with numerical simulations for models of multiclass traffic flow and polydisperse sedimentation.