2

I have a set of points defined by their longitudes and latitudes in WGS84 that define a polygon that I need to calculate the area of.

I originally tried to approach this using pyproj and Shapely as follows:

import numpy as np
from shapely.geometry import shape
import pyproj

points = np.array([[-14.51024374, 58.52358359], [-14.46760657, 58.36331163], [-14.42947386, 58.3302364 ], [-14.38354315, 58.29630198],[-14.28953719, 58.33202542], [-14.25769178, 58.35584472], [-14.23289391, 58.42266725], [-14.20653403, 58.76555668], [-14.26124968, 58.85410627], [-14.35780894, 58.92380724], [-14.43655603, 58.88558921], [-14.50306931, 58.75881054]])
min_lat = min(points[:,1])
max_lat = max(points[:,1])
mid_lat = (min_lat + max_lat) / 2
mid_lon = (min(points[:,0]) + max(points[:,0]))/2

a = pyproj.Proj(f"+proj=aea +lat_1={min_lat} +lat_2={max_lat} +lat_3={mid_lat} +lon_0={mid_lon} +ellps=WGS84 +datum=WGS84")
x, y = a(points[:,0], points[:,1])
cop = {"type": "Polygon", "coordinates": [zip(x, y)]}
shape(cop).area

which returns 960,473,396 m^2. Plainmeter agrees with this, and I have imported the polygon into QGIS and used the Field Calculator to calculate the area, which also agrees.

The polygon as defined in the .kml is:

<Polygon>
    <tessellate>1</tessellate>
    <gx:altitudeMode>clampToGround</gx:altitudeMode>
    <outerBoundaryIs>
        <LinearRing>
            <coordinates>
            -14.51024374,58.52358359,0 -14.46760657,58.36331163,0 -14.42947386,58.3302364,0 -14.38354315,58.29630198,0 -14.28953719,58.33202542,0 -14.25769178,58.35584472,0 -14.23289391,58.42266725,0 -14.20653403,58.76555668,0 -14.26124968,58.85410627,0 -14.35780894,58.92380724,0 -14.43655603,58.88558921,0 -14.50306931,58.75881054,0          
            </coordinates>
        </LinearRing>
    </outerBoundaryIs>
</Polygon>

Google Earth Pro measures this area as 908,350,894 m^2.

Is there a reason for this, or have I missed something?

1

I believe I have now resolved this.

The polygon in the .kml file is missing a closing point (i.e. the same as the first point). Adding the first point at the end of the coordinates tag, importing to Google Earth Pro and measuring the area produces 955km^2, much closer to that produced by the python code.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.