Here is my guess as to the problem. In your first integral with respect to r you don't indicate that you included the Jacobian that will be necessary in conversion to polar coordinates. In this case it is just r. That is, in differential notation dx*dy = r*dr*dphi. You want the double integral of:
r*F(r,phi)
with respect to r and phi within the limits you have described for the annulus. With the units I see in your picture that would make an enormous difference.
It was not entirely clear how you handled the meshgrid values x and y as a check. Perhaps you could explain this more clearly if the above is not the problem. You are bound to have some error in this method because of the numerous little squares that fall partly inside and partly outside the annulus. You haven't stated how large the discrepancy is between your two methods.
