  Many methods for modelling spatial processes assume global smoothness
properties; such assumptions are often violated in practice. We introduce a
method for modelling spatial processes that display heterogeneity or contain
discontinuities. The problem of non-stationarity is dealt with by using a
combination of Voronoi tessellation to partition the input space, and a
separate Gaussian process to model the data on each region of the partitioned
space. Our method is highly flexible because we allow the Voronoi cells to form
relationships with each other, which can enable non-convex and disconnected
regions to be considered. In such problems, identifying the borders between
regions is often of great importance and we propose an adaptive sampling method
to gain extra information along such borders. The method is illustrated with
simulation studies and application to real data.
