Spatial interpolation is a method that can be used to estimate the unknown value (or parameter) at certain time or location. Besides of that, shape preserving interpolation especially positivity preserving interpolation play an important role for data visualization where the resulting interpolating surface must be positive everywhere. In this study, we use a triangulation based method consists of positivity preserving cubic Bézier triangular patches to estimate the unknown value of rainfall amount at spatial localization. This scheme is useful for spatial data interpolation compared to the mesh free methods which are required the optimization technique in order to ensure the positivity of basis functions. The data is triangulated by using Delaunay's triangulation algorithm and from there the triangular surfaces are constructed piecewisely comprising three local schemes that are blended together via convex combination to form a scattered interpolating surface. To preserve the positivity of the surface, the appropriate conditions for the Bézier ordinates are implemented. Validation of the propose schemes are conducted by using wellknown test function used in previous researches. Then, the best the scheme is selected to visualize a rainfall data distribution based on secondary data collected at various meteorology stations in Malaysia obtain from Malaysian Meteorology Department. Finally, we estimate the unknown rainfall amount at certain locations that lies inside the convex hull formed by all given stations. The finding shows that positivity preserving cubic triangular patches is able to preserve the positivity of rainfall distribution surface and has a significant application in estimating the rainfall amount at unallocated meteorology stations in Malaysia. © 2018, Insight Society.