经纬度轮廓面积计算

import js2py
import pandas as pd
import numpy as np
from shapely.geometry import Polygon,Point


data_js = open(r'polygon_go.js', 'r', encoding='utf8').read()
data_js = js2py.eval_js(data_js)

def def_get_area_js(poly,data_js):
    '''
    计算poly的面积(单位：平方公里)
    :param :poly是一个polygon或者是一个字符串（wkt）
    :return :返回polygon的面积,若无polygon,默认为0
    '''
    try:
        shape = str(poly).replace(', ',',')
        shape = shape.replace('POLYGON ((','').replace('))','').replace(', ',',')
        k2 = shape.split(',')
        k3 = []
        for i in k2:
            print(i)
            i = i.replace(")","").replace("(","")
            lng_1 = float(i.split(' ')[0])
            lat_1 = float(i.split(' ')[1])
            k3.append([lng_1,lat_1])

        final = data_js(k3)
        poly_area = float(final)*0.000001
        return poly_area
    except:
#        print()
        return 0
   
data = pd.read_excel(')

id_area = []
for index,row in data.iterrows():
    # wkt = '108.3886230000002 22.73205900000006, 108.3749951436175 22.73057592190073...'
    wkt = row['轮廓坐标串']
    wkt = str(wkt).replace('POLYGON ((','').replace('))','').replace(', ',',').replace('\n',',')## 一般不需要用到
    shape = wkt.split(',')
    shape = np.array([i.split(' ') for i in shape])
    shape = shape.astype(np.float)
    poly = Polygon(shape)
    if poly.is_valid == True:
        pass
    else:
        poly = poly.buffer(0)

    area = def_get_area_js(poly,data_js)

用到的 js 文件代码

    function getArea(ring) {
        var sJ = 6378137;
        var Hq = 0.017453292519943295;

        var arr = [];
        for (var i = 0;i< ring.length; i++){
            var poi = {
                lng:ring[i][0],
                lat:ring[i][1]
            };
            arr[i] = poi;
        }

        ring = arr;

        var c = sJ *Hq , d = 0 , e = ring.length;

        if (3 > e) {
            return 0;
        }

        for (var g = 0; g < e - 1; g += 1){
            var h = ring[g], k = ring[g + 1];
            var u = h.lng * c * Math.cos(h.lat * Hq);
            var h = h.lat * c;
            var v = k.lng * c * Math.cos(k.lat *Hq);
            var d = d + (u * k.lat * c - v * h);
        }

        g = ring[g];
        ring = ring[0];
        e = g.lng * c * Math.cos(g.lat * Hq);
        g = g.lat * c;
        k = ring.lng * c * Math.cos(ring.lat * Hq);
        d += e * ring.lat * c - k * g;
        return 0.5*Math.abs(d)
    }