Accompanying source code for my article: https://erkaman.github.io/posts/area_convex_polygon.html

shoelace.cpp

/*
This software is released under the MIT license:

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/

#include <stdio.h>

/*
  It is assumed that the vertex coordinates form a convex polygon,
  and that the vertices are given in a counter-clockwise order.
  We do no error checking, to keep things simple.
*/
float area(float *x, float *y, int m) {
    int n = m-1;
    float s = 0.0f;
    for(int i = 0; i <= n; ++i) {
        s += x[i] * y[(i+1)%(n+1)] - x[(i+1)%(n+1)] * y[i];
    }
    return 0.5f * s;
}

int main() {
    /*
      We use the derived shoelace formula to compute some areas.
     */

    {
        const int N = 3;
        float x[N] = {+9.0f, +7.0f, +5.0f};
        float y[N] = {+4.0f, +8.0f, +4.0f};
        // hand calculation gives that the expected area is 8.0
        printf("area of triangle %f\n", area(x, y, N));
    }

    {
        const int N = 4;
        float x[N] = {+4.0f, -4.0f, -4.0f, +4.0};
        float y[N] = {+2.0f, +2.0f, -2.0f, -2.0};

        // hand calculation gives that the expected area is 32.0
        printf("area of rectangle %f\n", area(x, y, N));
    }

    {
        const int N = 5;
        float x[N] = {+13.0f, +9.0f, +5.0f, +7.0f, +11.0f};
        float y[N] = {-3.0f, +1.0f, -3.0f, -7.0f, -7.0f};

        // hand calculation gives that the expected area is 40.0
        printf("area of pentagon %f\n", area(x, y, N));
    }
}

I believe there's a numerical instability in the formula: if a polygon is quite far away from coordinates origin, but quite small in size, x_i * y _j - x_j * y_i can suffer from catastrophic cancellation and loose precision. In my project I've tried to overcome that by subtracting minimal X and Y values from corresponding coordinates (written is TypeScript):

/**
 * Computes signed area of a ngon.
 * @see http://mathworld.wolfram.com/PolygonArea.html
 *
 * @param ngon Set of points.
 * @param start Index of the first point of the ngon in the given set.
 * @param end Index of the last point of the ngon in the given set.
 * @returns Signed area. For simple polygons it's less than zero for clockwise
 *      ordered ngons and greater than zero for counterclockwise ones.
 */
export default function signedArea(
    ngon: ArrayLike<number>,
    start: number = 0,
    end: number = ngon.length
): number {
    // There's possibility of losing precision here due to multiplication of large
    // numbers. "Shifting" the whole ngon by coordinates of lower-left corner of
    // its bbox somewhat alleviate the problem, but not for large enough ngons.
    // Still, we can check winding order.

    let minX = ngon[start];
    let minY = ngon[start + 1];

    for (let i = start + 2; i < end; i += 2) {
        minX = ngon[i] < minX ? ngon[i] : minX;
        minY = ngon[i + 1] < minY ? ngon[i + 1] : minY;
    }

    let area = 0;

    for (let i = end - 2, j = start; j < end; i = j, j += 2) {
        area += ((ngon[i] - minX) * (ngon[j + 1] - minY) -
            (ngon[j] - minX) * (ngon[i + 1] - minY)) / 2;
    }

    return area;
}

@dmikis - Why not https://en.wikipedia.org/wiki/Trapezoidal_rule? Works for any polygon.

function getSquare (points) {
    let ret = 0;
    for (const i = 0, l = countor.length - 1; i < l; ++i) {
        ret += (points[i][0] - points[i + 1][0]) * (points[i][1] + points[i + 1][1] /* - minY ? */) / 2;
    }
    return ret;
}