Validating and controlling safety-critical systems in uncertain environments necessitates probabilistic reachable sets of future state evolutions. The existing methods of computing probabilistic reachable sets normally assume that stochastic uncertainties are independent of system states, inputs, and other environment variables. However, this assumption falls short in many real-world applications, where the probability distribution governing uncertainties depends on these variables, referred to as contextual uncertainties. This paper addresses the challenge of computing probabilistic reachable sets of stochastic nonlinear states with contextual uncertainties by seeking minimum-volume polynomial sublevel sets with contextual chance constraints. The formulated problem cannot be solved by the existing sample-based approximation method since the existing methods do not consider conditional probability densities. To address this, we propose a consistent sample approximation of the original problem by leveraging conditional density estimation and resampling. The obtained approximate problem is a tractable optimization problem. Additionally, we prove the proposed sample-based approximation’s almost uniform convergence, showing that it gives the optimal solution almost consistently with the original ones. Through a numerical example, we evaluate the effectiveness of the proposed method against existing approaches, highlighting its capability to significantly reduce the bias inherent in sample-based approximation without considering a conditional probability density.