Immediately after a major earthquake strikes, a segment of the people in the disaster region must evacuate to relief shelters. If the evacuation proceeds in an efficient and safe manner, many lives can potentially be saved as the effects of secondary factors, e.g., secondary collapses, stampedes, and aftershocks, can largely be mitigated, while at the same time the people taking refuge can quickly gain access to vital relief goods. Formulating the best possible protocols for post-earthquake evacuation is a crucial, yet highly challenging problem considering the profound uncertainty which is present in the post-disaster environment, including the location and degree of infrastructural damage, the number and location of people who need to evacuate, the behavior and psychology of evacuees, etc. In this research, we formulate a two-stage stochastic evacuation simulation model capable of credibly accounting for these uncertainties. We develop a solution methodology to solve the shelter location–allocation problem represented by the two-stage model, finding the optimal location of shelters to be open as the first stage solution and the optimal assignment of evacuees to shelters as the second stage solution to minimize the total expected cost. The cost-based objective function consists of the cost of opening shelters, a deprivation cost component, the cost of trapped evacuees, and the cost of shelters exceeding capacity. In developing this framework, we collaborate with the National Science and Technology Center for Disaster Reduction (NCDR) in utilizing GIS-based probabilistic models of road damage. We apply the proposed framework to run an extensive numerical study and present various valuable observations and insights which can be gained.