The analytic solution of the gravity anomaly caused by a 2D irregular mass body with the density contrast varying as a polynomial function in the horizontal and vertical directions is extrapolated from a historical version in which the analytic solution for the gravity anomaly was given only at the origin of the coordinate system to any point for the density function in terms of variables relative to that origin. To calculate the gravity anomaly at stations that are not at origins, a coordinate transformation is performed, in which case the polynomial density contrast function must also be expressed in the transformed coordinates, or a transformed solution must be obtained. These analytic solutions can be obtained at any station using (1) a solution transformation method, in which the density function and boundary of a mass body are kept intact, or (2) a coordinate transformation method, in which polynomial coefficient and boundary of a mass body are transformed accordingly. The issue of singularity and instability of the analytic methods has been related to case studies. Caution should be exercised in modeling or interpreting the gravity survey data using the analytic methods for large target-distance-to-target-size ratios outside the range of numerical stability. Compared with other published methods, the analytic solution results agree very well with other numerical or seminumerical methods, indicating the solution is correct and can be applied for any gravity anomaly calculation caused by an irregular 2D mass body with the density-contrast approximated as a polynomial function of horizontal position and/or vertical position when the observation is within the range of numerical stability.