Climatic response of population boundaries, known as ecotone, is a classic indicator in ecology. In addition to known challenges, the spatial and dynamical characteristics of the boundary are not only affected by the spatial gradient in the environmental factors, but also by local heterogeneities in the regional characteristics. Here, we capture the effects of time-independent, quenched spatial heterogeneities on the ecological boundary with the disordered contact process in one and two dimensions with a linear spatial trend in the local control parameter. We apply the strongdisorder renormalization group method to calculate the sites occupied with an O(1) probability in the stationary state, readily yielding the population front's position as the outermost site locally as well as globally for the entire boundary. We show that under a quasistatic change of the global environment, mimicking climate change, the front advances intermittently: long quiescent periods are interrupted by rare but long jumps. The characteristics of this intermittent dynamics are found to obey universal scaling laws in terms of the gradient, conjectured to be related to the correlation-length exponent of the model. In the two dimensional case our analysis shows that both the local and global shift of the front exhibit intermittent dynamics and follow the same scaling laws as in one dimension, apart from a logarithmic correction for global shifts. In one dimension, we show that determining the front position is related to an extremal problem of random walks, shedding further light on the origin of the jumps. Our results are in stark contrast to the behavior of the model without disorder, suggesting that current observations might misleadingly show little to no climate response for an extended period of time, concealing the long-term effects of climate change.