2014-02-26 184 views
7

我正在嘗試爲MKPolygon製作面積計算類別。 我發現了一些JS代碼https://github.com/mapbox/geojson-area/blob/master/index.js#L1,並鏈接到算法:http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409。 它說: enter image description hereMKPolygon面積計算

enter image description here

這裏是我的代碼,這給了一個錯誤的結果(比實際數千倍以上):

#define kEarthRadius 6378137 
@implementation MKPolygon (AreaCalculation) 
- (double) area { 
    double area = 0; 
    NSArray *coords = [self coordinates]; 
    if (coords.count > 2) { 
     CLLocationCoordinate2D p1, p2; 
     for (int i = 0; i < coords.count - 1; i++) { 
      p1 = [coords[i] MKCoordinateValue]; 
      p2 = [coords[i + 1] MKCoordinateValue]; 
      area += degreesToRadians(p2.longitude - p1.longitude) * (2 + sinf(degreesToRadians(p1.latitude)) + sinf(degreesToRadians(p2.latitude))); 
     } 

     area = area * kEarthRadius * kEarthRadius/2; 
    } 

    return area; 
} 
- (NSArray *)coordinates { 
    NSMutableArray *points = [NSMutableArray arrayWithCapacity:self.pointCount]; 
    for (int i = 0; i < self.pointCount; i++) { 
     MKMapPoint *point = &self.points[i]; 
     [points addObject:[NSValue valueWithMKCoordinate:MKCoordinateForMapPoint(* point)]]; 
    } 
    return points.copy; 
} 

double degreesToRadians(double radius) { 
    return radius * M_PI/180; 
} 

@end 

我錯過了什麼?

+0

您的循環中缺少'i == N-1'和'i + 1 == 0'(環繞)的最後一步。 –

+0

@MartinR謝謝,我會解決它。 – Shmidt

+0

@MartinR你是對的。如果您發表評論作爲答案,我會接受它。 – Shmidt

回答

3

循環中缺少i = N-1i+1 = 0(環繞)的最後一步。

+0

有人能解釋這實際上代碼形式的含義嗎? – AndyDunn

1

這可以幫助別人...... 你需要的形狀邊緣點傳遞到下面的方法,它返回一個多邊形雨燕3.0實施

static double areaOfCurveWithPoints(const NSArray *shapeEdgePoints) { 

    CGPoint initialPoint = [shapeEdgePoints.firstObject CGPointValue]; 

    CGMutablePathRef cgPath = CGPathCreateMutable(); 
    CGPathMoveToPoint(cgPath, &CGAffineTransformIdentity, initialPoint.x, initialPoint.y); 

    for (int i = 1;i<shapeEdgePoints.count ;i++) { 

     CGPoint point = [[shapeEdgePoints objectAtIndex:i] CGPointValue]; 
     CGPathAddLineToPoint(cgPath, &CGAffineTransformIdentity, point.x, point.y); 
    } 
    CGPathCloseSubpath(cgPath); 

    CGRect frame = integralFrameForPath(cgPath); 
    size_t bytesPerRow = bytesPerRowForWidth(frame.size.width); 
    CGContextRef gc = createBitmapContextWithFrame(frame, bytesPerRow); 
    CGContextSetFillColorWithColor(gc, [UIColor whiteColor].CGColor); 
    CGContextAddPath(gc, cgPath); 
    CGContextFillPath(gc); 

    double area = areaFilledInBitmapContext(gc); 

    CGPathRelease(cgPath); 
    CGContextRelease(gc); 

    return area; 
} 
static CGRect integralFrameForPath(CGPathRef path) { 
    CGRect frame = CGPathGetBoundingBox(path); 
    return CGRectIntegral(frame); 
} 

static size_t bytesPerRowForWidth(CGFloat width) { 
    static const size_t kFactor = 64; 
    // Round up to a multiple of kFactor, which must be a power of 2. 
    return ((size_t)width + (kFactor - 1)) & ~(kFactor - 1); 
} 

static CGContextRef createBitmapContextWithFrame(CGRect frame, size_t bytesPerRow) { 
    CGColorSpaceRef grayscale = CGColorSpaceCreateDeviceGray(); 
    CGContextRef gc = CGBitmapContextCreate(NULL, frame.size.width, frame.size.height, 8, bytesPerRow, grayscale, kCGImageAlphaNone); 
    CGColorSpaceRelease(grayscale); 
    CGContextTranslateCTM(gc, -frame.origin.x, -frame.origin.x); 
    return gc; 
} 

static double areaFilledInBitmapContext(CGContextRef gc) { 
    size_t width = CGBitmapContextGetWidth(gc); 
    size_t height = CGBitmapContextGetHeight(gc); 
    size_t stride = CGBitmapContextGetBytesPerRow(gc); 

    // Get a pointer to the data 
    unsigned char *bitmapData = (unsigned char *)CGBitmapContextGetData(gc); 

    uint64_t coverage = 0; 
    for (size_t y = 0; y < height; ++y) { 
     for (size_t x = 0; x < width; ++x) { 
      coverage += bitmapData[y * stride + x]; 
     } 
    } 
    // NSLog(@"coverage =%llu UINT8_MAX =%d",coverage,UINT8_MAX); 
    return (double)coverage/UINT8_MAX; 
} 
2

整個算法的正確區域:

import MapKit 
let kEarthRadius = 6378137.0 

// CLLocationCoordinate2D uses degrees but we need radians 
func radians(degrees: Double) -> Double { 
    return degrees * M_PI/180; 
} 

func regionArea(locations: [CLLocationCoordinate2D]) -> Double { 

    guard locations.count > 2 else { return 0 } 
    var area = 0.0 

    for i in 0..<locations.count { 
     let p1 = locations[i > 0 ? i - 1 : locations.count - 1] 
     let p2 = locations[i] 

     area += radians(degrees: p2.longitude - p1.longitude) * (2 + sin(radians(degrees: p1.latitude)) + sin(radians(degrees: p2.latitude))) 
    } 

    area = -(area * kEarthRadius * kEarthRadius/2); 

    return max(area, -area) // In order not to worry about is polygon clockwise or counterclockwise defined. 
}