Math Behind MKMetersPerMapPointAtLatitude

Question

I'm trying to convert some Apple mapping code to Java. I have most of it converted correctly except for a few calls to MKMetersPerMapPointAtLatitude

I have a very close solution... but it's not exact and I'm not sure why not. Any ideas?

#import <Foundation/Foundation.h>
#import <Math.h>
@import MapKit;

#define MERCATOR_OFFSET 268435456.0 / 2.0
#define MERCATOR_RADIUS (MERCATOR_OFFSET/M_PI)
#define WGS84_RADIUS 6378137.0
#define POINTS_PER_METER (MERCATOR_RADIUS / WGS84_RADIUS)

double MyMetersPerMapPointAtLatitude(double latitude) {
    return 1.0 / (POINTS_PER_METER / cos(latitude * M_PI / 180.0));
}

int main(int argc, const char * argv[]) {
    @autoreleasepool {
        double latitude = 33.861315;
        for (int i = 0; i < 100; i++) {
            double a = MKMetersPerMapPointAtLatitude(latitude);
            double b = MyMetersPerMapPointAtLatitude(latitude);

            NSLog(@"%f %f", a, b);
            latitude += .1;
        }
    }
    return 0;
}

Prints out the following

2015-05-19 09:13:00.334 Test[92619:5369062] 0.123522 0.123969
2015-05-19 09:13:00.335 Test[92619:5369062] 0.123379 0.123824
2015-05-19 09:13:00.335 Test[92619:5369062] 0.123236 0.123678
2015-05-19 09:13:00.335 Test[92619:5369062] 0.123092 0.123532
2015-05-19 09:13:00.335 Test[92619:5369062] 0.122948 0.123386
2015-05-19 09:13:00.335 Test[92619:5369062] 0.122804 0.123239
2015-05-19 09:13:00.335 Test[92619:5369062] 0.122659 0.123092
...etc

Randolph

What exactly you think is not working? – Tate

Well it's not quite accurate but I'm not sure why. It appears my function is nearly accurate if I add 0.3 to the latitude... – Randolph

At a guess, something to do with the fact that the Earth is an oblate ellipsoid, not a perfect sphere? It's fatter around the equator. – Impurity

Answer 1

For starters, we can rearrange your function to make it a bit more readable:

double MyMetersPerMapPointAtLatitude(double latitude) {
    return cos(latitude * M_PI / 180.0) / POINTS_PER_METER;
}

Now, the issue is, as Tommy points out, that you're not accounting for the flattening of the Earth. You can do that with:

double f = 1/298.257223563; // WGS84 flattening
double MyMetersPerMapPointAtLatitude(double latitude) {
    return (1-f) * cos(latitude * M_PI / 180.0) / POINTS_PER_METER;
}

That gets the error down to the point where I would attribute it to rounding and truncation.

Forequarter