Calculating Geometric Offsets Based on Target Area in GIS

Posted on Edited on In

Introduction

Assume that all coordinates are expressed in a suitable projected Cartesian coordinate system, so Euclidean lengths and planar areas are meaningful. Land subdivision sometimes involves moving a line segment between two non-parallel boundaries to produce a polygon with a specified area.

The Scenario:
A user initiates a tool and defines:

  1. Line 1: A fixed vertex A and a direction point B.
  2. Line 2: A fixed vertex C and a direction point D.
  3. Target Area: The required area of the generated trapezoid, expressed in the square units of the projected coordinate system.

The Goal:
Calculate the offset distance and the coordinates of the new segment such that the area of the polygon equals the Target Area.

The Geometric Concept: The “Virtual Triangle”

Because lines and are not parallel, their infinite extensions intersect at a unique point. Let us call this intersection point .

  1. The Base Triangle: The segments and form a triangle (the “Virtual Triangle”).
  2. The Scaled Triangle: Moving to creates another triangle, , which may be larger or smaller than .
  3. Similarity: Multiplying both vectors from by the same scale factor makes parallel to . Therefore, is similar to .

The area of the desired trapezoid () is the absolute difference between the areas of the two triangles.

The Scaling Law

For similar triangles, the ratio of their areas is equal to the square of the ratio of their corresponding dimensions (side lengths or heights).

If we let be the area of and be the area of , the non-negative scaling factor is:

Once is found, the new coordinates and can be found by simply scaling the vectors starting from the intersection point .

The Derivation

To solve this programmatically, we need to derive the formula for the offset distance .

Let:

  • be the area of the virtual triangle .
  • be the length of segment .
  • be the altitude of from to base .
  • be the user-specified target area.

From triangle geometry, we know:

For an outward subdivision, and the new height is . The scaling law gives

Solving for :

Substituting :

For an inward subdivision, and the new height is . Therefore,

The Algorithm

To implement this in a GIS tool, such as a QGIS plugin or an ArcPy script, follow this logic:

  1. Find Intersection ():
    Compute the intersection point of the infinite lines defined by vectors and . Use a scale-aware tolerance to identify parallel or nearly parallel lines.
    Parallel boundaries require a separate construction and are outside the non-parallel scenario considered here.

  2. Calculate Initial Virtual Area ():
    Calculate the area of the triangle formed by .

  3. Determine Direction (Add or Subtract Area):
    We must determine if the user wants to expand the triangle (offset away from ) or shrink it (offset toward ).

    • Calculate both and .
    • Both dot products must have the same sign; otherwise, the two direction points select inconsistent branches from .
    • If both are positive: The boundaries point away from , so .
    • If both are negative: The boundaries point toward , so .

  4. Calculate Scale Factor ():

  5. Compute New Coordinates:
    Scale the vectors from the origin :

  6. Compute Height ():
    Calculate the perpendicular distance from to the infinite line through and .

Python Implementation

Below is a simple Python function that performs this calculation.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import math


def calculate_offset_geometry(
    ptA, ptB, ptC, ptD, target_area, epsilon=1e-8
):
    """
    Calculate offset distance h and new coordinates A', C' from a target area.

    All points must use the same projected Cartesian coordinate system.

    Args:
        ptA, ptB: (x, y) tuples defining line 1.
        ptC, ptD: (x, y) tuples defining line 2.
        target_area: Required area of polygon AA'C'C.
        epsilon: Relative tolerance for geometric comparisons.

    Returns:
        A dictionary containing h, A_prime, and C_prime.
    """
    if epsilon <= 0:
        raise ValueError("epsilon must be positive.")
    if not math.isfinite(target_area) or target_area < 0:
        raise ValueError("target_area must be a finite, non-negative number.")

    coordinates = (*ptA, *ptB, *ptC, *ptD)
    if len(coordinates) != 8 or not all(map(math.isfinite, coordinates)):
        raise ValueError("Each point must contain two finite coordinates.")

    def vector(start, end):
        return (end[0] - start[0], end[1] - start[1])

    def dot(u, v):
        return u[0] * v[0] + u[1] * v[1]

    def cross(u, v):
        return u[0] * v[1] - u[1] * v[0]

    def norm(v):
        return math.hypot(v[0], v[1])

    direction_ab = vector(ptA, ptB)
    direction_cd = vector(ptC, ptD)
    base = vector(ptA, ptC)
    length_ab = norm(direction_ab)
    length_cd = norm(direction_cd)
    base_length = norm(base)
    geometry_scale = max(1.0, length_ab, length_cd, base_length)
    length_tolerance = epsilon * geometry_scale

    if length_ab <= length_tolerance or length_cd <= length_tolerance:
        raise ValueError("Each boundary line requires two distinct points.")
    if base_length <= length_tolerance:
        raise ValueError("A and C must be distinct points.")

    denominator = cross(direction_ab, direction_cd)
    if abs(denominator) <= epsilon * length_ab * length_cd:
        raise ValueError("Boundary lines are parallel or nearly parallel.")

    # A + parameter * AB is the intersection O.
    parameter = cross(vector(ptA, ptC), direction_cd) / denominator
    O = (
        ptA[0] + parameter * direction_ab[0],
        ptA[1] + parameter * direction_ab[1],
    )

    vector_oa = vector(O, ptA)
    vector_oc = vector(O, ptC)
    length_oa = norm(vector_oa)
    length_oc = norm(vector_oc)
    twice_initial_area = abs(cross(vector_oa, vector_oc))

    if (
        length_oa <= length_tolerance
        or length_oc <= length_tolerance
        or twice_initial_area <= epsilon * length_oa * length_oc
    ):
        raise ValueError("The virtual triangle is degenerate.")

    direction_a = dot(vector_oa, direction_ab)
    direction_c = dot(vector_oc, direction_cd)
    ambiguous_a = epsilon * length_oa * length_ab
    ambiguous_c = epsilon * length_oc * length_cd

    if abs(direction_a) <= ambiguous_a or abs(direction_c) <= ambiguous_c:
        raise ValueError("A boundary direction is numerically ambiguous.")
    if (direction_a > 0) != (direction_c > 0):
        raise ValueError("AB and CD select inconsistent branches from O.")

    initial_area = twice_initial_area / 2.0
    if direction_a > 0:
        new_area = initial_area + target_area
    else:
        if target_area > initial_area and not math.isclose(
            target_area, initial_area, rel_tol=epsilon
        ):
            raise ValueError(
                "Target area exceeds the size of the converging triangle tip."
            )
        new_area = max(0.0, initial_area - target_area)

    scale = math.sqrt(new_area / initial_area)
    A_prime = (
        O[0] + vector_oa[0] * scale,
        O[1] + vector_oa[1] * scale,
    )
    C_prime = (
        O[0] + vector_oc[0] * scale,
        O[1] + vector_oc[1] * scale,
    )

    displacement = vector(ptA, A_prime)
    h = abs(cross(base, displacement)) / base_length

    return {
        "h": h,
        "A_prime": A_prime,
        "C_prime": C_prime,
    }

Another Perspective: Determinants and Vector Projection

The virtual-triangle method can also be described entirely using linear algebra. This perspective treats triangle area as a determinant and the offset distance as the length of a vector component. It provides a compact alternative formulation and another way to verify the geometry. In this section, is the same scale factor as above, and is the same target area.

Translate the coordinate system so that the intersection point is the origin. Define

The determinant of is the signed area of the parallelogram formed by and . Therefore, the area of the virtual triangle is

Scale both vectors by the same factor :

The segment joining the scaled endpoints is therefore parallel to .

If , then

For an outward subdivision with target area ,

Substituting gives

For an inward subdivision, the new triangle is smaller, so

Finally, translate the scaled vectors back to the original coordinate system:

Original vectors and their scaled versions forming a parallel subdivision

Finding the offset distance by projection

Let the displacement from to be

and let the original base direction be

The component of parallel to is

Subtracting this projection leaves the perpendicular component. Its norm is the offset distance:

NumPy implementation

The following function applies this formulation after the intersection point has been found. Passing outward=False selects the inward solution.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
import math
import numpy as np


def subdivision_by_linear_algebra(
    ptO, ptA, ptC, target_area, outward=True, epsilon=1e-12
):
    if epsilon <= 0:
        raise ValueError("epsilon must be positive.")
    if not math.isfinite(target_area) or target_area < 0:
        raise ValueError("target_area must be a finite, non-negative number.")

    O = np.asarray(ptO, dtype=float)
    A = np.asarray(ptA, dtype=float)
    C = np.asarray(ptC, dtype=float)

    points = (O, A, C)
    if any(point.shape != (2,) for point in points):
        raise ValueError("Each point must contain exactly two coordinates.")
    if not all(np.all(np.isfinite(point)) for point in points):
        raise ValueError("All coordinates must be finite.")

    a = A - O
    c = C - O
    base = C - A
    length_a = np.linalg.norm(a)
    length_c = np.linalg.norm(c)
    base_length = np.linalg.norm(base)
    geometry_scale = max(1.0, length_a, length_c, base_length)
    length_tolerance = epsilon * geometry_scale
    matrix = np.column_stack((a, c))
    twice_initial_area = abs(np.linalg.det(matrix))

    if (
        length_a <= length_tolerance
        or length_c <= length_tolerance
        or twice_initial_area <= epsilon * length_a * length_c
    ):
        raise ValueError("The virtual triangle is degenerate.")

    sign = 1.0 if outward else -1.0
    scale_squared = 1.0 + sign * 2.0 * target_area / twice_initial_area

    if scale_squared < -epsilon:
        raise ValueError("The inward target area exceeds the virtual triangle.")
    scale = math.sqrt(max(0.0, scale_squared))

    A_prime = O + scale * a
    C_prime = O + scale * c

    displacement = A_prime - A
    base_squared = np.dot(base, base)
    if base_squared <= length_tolerance**2:
        raise ValueError("A and C must be distinct points.")

    parallel_component = np.dot(displacement, base) / base_squared * base
    perpendicular_component = displacement - parallel_component
    h = np.linalg.norm(perpendicular_component)

    scaled_matrix = scale * matrix
    calculated_area = abs(
        abs(np.linalg.det(scaled_matrix)) - twice_initial_area
    ) / 2.0

    return {
        "scale": scale,
        "h": float(h),
        "A_prime": tuple(map(float, A_prime)),
        "C_prime": tuple(map(float, C_prime)),
        "calculated_area": float(calculated_area),
    }

For example, take , , , and :

result = subdivision_by_linear_algebra(
    ptO=(0, 0),
    ptA=(8, 1),
    ptC=(6, 3),
    target_area=20,
)

print(result)

This produces approximately:

{
    'scale': 1.79505494,
    'h': 5.05969863,
    'A_prime': (14.36043949, 1.79505494),
    'C_prime': (10.77032961, 5.38516481),
    'calculated_area': 20.0
}

Here, and . Their difference is , so the area between the two parallel segments is , as required.