Mapping of disease incidence has long been of importance to epidemiology and public health. In this paper, we consider identification of clusters of spatial units with elevated disease rates and develop a new approach that estimates the relative disease risk in association with potential risk factors and simultaneously identifies clusters corresponding to elevated risks. A heterogeneity measure is proposed to enable the comparison of a candidate cluster and its complement under a pair of complementary models. A quasi-likelihood procedure is developed for estimating the model parameters and identifying the clusters. An advantage of our approach over traditional spatial clustering methods is the identification of clusters that can have arbitrary shapes due to abrupt or non-contiguous changes while accounting for risk factors and spatial correlation. Asymptotic properties of the proposed methodology are established and a simulation study shows empirically sound finite-sample properties. The mapping and clustering of enterovirus 71 infection in Taiwan are carried out for illustration.