The probability integral method is an important approach to predicting mining subsidence, which can be used to calculate surface movement and deformation with normal distribution profile in the gently inclined and inclined coal seams. To overcome the shortcoming of the original method for arbitrary shape face, the integral regions of mining are defined by the strike and dip directions with different integral functions, and divided into several unstructured triangular elements as the basic unit of double integrals. By the double integral conversion method with coordinate rotation transformation, the integral of the basic unit can be converted into double integrals with two lines for the upper and lower limits. When the influence radius of rotation is given, the new double integrals can be calculated using the numerical compound Simpson method to obtain the subsidence at any point on the surface due to the mining of basic triangular element. Finally, the superposition calculation of all basic triangular elements is used to finish the probability integral at any point to predict mining subsidence for arbitrary shape face. The proposed algorithm is implemented in the GIS platform with good results in the case study, providing the support for “under three” coal mining with surface movement prediction.