The recent advances in remote sensing technology resulted in peta bytes of data in raster format. To process this data, it is often combined with high resolution vector data that represents, for example, region boundaries. One of the common operations that combine big vector and raster data is the zonal statistics which computes some aggregate values for each polygon in the vector dataset. This paper proposes a novel and scalable algorithm for zonal statistics that can scale to peta bytes of raster and vector data. The proposed method does not require any preprocessing or indexing making it perfect for ad-hoc queries that scientists usually want to run. We implement a prototype for the proposed method and the initial preliminary results show that the proposed method can scale up-to a trillion pixels.