“Is this point to the left or to the right of my path?” Or a more intuitive situation: Is the address point to the left or right of the street? This post explores how to solve this using the 3D cross product with NumPy and Shapely libraries.

1. The Core Problem: Side Determination

When moving along a directed line from point to point , any point in the 2D plane must fall into one of three categories:

1. Left Side: A counter-clockwise turn is needed to face the point.
2. Right Side: A clockwise turn is needed to face the point.
3. Collinear: The point lies exactly on the line (or its infinite extension).

To solve this mathematically, we treat the segments as vectors. Let be our direction and be the offset to our target point.

2. The Math: The 3D Cross Product

While the cross product is natively a 3D operation, it is the perfect tool for 2D side determination.

Geometric Definition

The 3D cross product () produces a new vector that is perpendicular to both inputs. By placing our 2D vectors on the plane (setting ), the resulting cross product vector must point directly along the -axis.

The Right-Hand Rule

The direction of this value follows the Right-Hand Rule:

• If your fingers curl from to and your thumb points Up (), the point is to the Left.
• If your thumb points Down (), the point is to the Right.

Algebraic Formula

For 2D vectors extended to 3D, the calculation simplifies significantly:

3. Applied Scenario: Path Following

In a real-world application, we don’t just have a single segment; we have a Polyline (). To determine the side relative to the “flow” of the path, we follow these steps:

1. Find the Nearest Point (): Locate the point on the polyline closest to our target (using Shapely’s project method).
2. Establish Direction (): * Normal Case: Walk forward along the polyline from to find point .
   • Endpoint Case: If is the end of the line, walk backward to find a point . To maintain a consistent forward-facing vector, we swap them: let the original be and the new be .
3. Construct Vectors: Define and .

4. Implementation: NumPy & Shapely

Using Shapely for linear referencing and NumPy for the vector math allows for concise, readable code. We use an epsilon threshold to account for floating-point inaccuracies.

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

import numpy as np
from shapely.geometry import LineString, Point

def get_side_of_path(L: LineString, A_geom: Point, distance: float = 5.0, epsilon: float = 1e-9):
    # 1. Find Point B (Nearest point on L)
    dist_to_B = L.project(A_geom)
    B_geom = L.interpolate(dist_to_B)
    
    line_length = L.length
    
    # 2. Logic for Point C (Look-ahead or Endpoint Flip)
    if dist_to_B < line_length:
        # Standard: Walk forward 5m (clamped to end)
        dist_to_C = min(dist_to_B + distance, line_length)
        C_geom = L.interpolate(dist_to_C)
        B_vec = np.array([B_geom.x, B_geom.y, 0])
        C_vec = np.array([C_geom.x, C_geom.y, 0])
    else:
        # Endpoint: Walk backward 5m, then flip B and C
        dist_to_C_prime = max(line_length - distance, 0)
        C_prime_geom = L.interpolate(dist_to_C_prime)
        C_vec = np.array([B_geom.x, B_geom.y, 0]) 
        B_vec = np.array([C_prime_geom.x, C_prime_geom.y, 0])

    # 3. Vector math
    A_vec = np.array([A_geom.x, A_geom.y, 0])
    u = C_vec - B_vec
    v = A_vec - B_vec
    
    # 4. Cross Product Calculation
    z_val = np.cross(u, v)[2]
    
    # 5. Interpretation with Epsilon
    if abs(z_val) < epsilon:
        return 0, "Collinear"
    elif z_val > 0:
        return 1, "Left"
    else:
        return -1, "Right"

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.

Question: Offset one edge of a quadrilateral (giving a new edge ) to form a new (usually smaller) quadrilateral (parcel) whose area equals a predefined value . Find the offset distance .

Short answer: find the linear equation between offset distance and offset edge , then solve the quadratic equation for the offset :

It is easy to pick the correct root as the solution for .

The linear relationship between and follows from trigonometry:

The quadratic equation of is:

There are other ways to explain .

Let be the origin; , and are vectors, so , , where is the scale ratio.
Then the linear equation relating and can be proved.

From the perspective of similar triangles,

In a GIS work environment, offsetting a line or reading the length of a line is straightforward and convenient. In order to get the equation of and , a trick can be used here.

Start from the general form:

Let ; then

Offset by a specific value — say 10 — and read from the GIS software; then

The problem is solved without needing angles or the distance , which could be far away.

My practice is to put the logic into an Excel form with three parameters.

Area of the new parcel, in
Length of the original edge, in
Length of the edge after offsetting by , in
Required offset, in B69: =IFERROR((-B67+SQRT(B67B67+2(B68-B67)/10*B66))/((B68-B67)/10),0)

The Openwrt VM running in my fanless pc is the main router of home network. Recently I got a new router which has 6 NICs and a more powerful CPU. So I need to migrate the Openwrt VM to the new machine.

The basic idea of migration is backup the vm, move the omv file to the new machine, then restore it. Here are the steps for the whole process.

1. Install Proxmox VE on the new machine and set up all the network bridges.
2. Backup the Openwrt VM, move the backup file to the new Proxmox VE node, and restore the VM by following the official wiki Backup_and_Restore.
3. On the new PVE node, remove the old network devices, then add new network devices according to the new hardware platform and Linux bridge setting.
4. Shut down the old Openwrt VM.
5. Start the new Openwrt VM, login to web ui. Update the device list for WAN and BR-LAN.
6. Connect the WAN cable to the designated NIC on the new PVE node, restart the VM and the modem if necessary. Then the new Openwrt VM should work good as router.

1. Introduction

In offshore petroleum exploration practice, the coordinates of blocks designed by geologists for a particular purpose must be provided in Degree-Minute-Second format: the degree uses 2 or 3 digits, the minute uses 2 digits, and the second uses 2 digits, in one of four values — 00, 15, 30 or 45 (in increments of a quarter of a minute). All coordinates must be listed in clockwise order, starting from the name vertex (the most northeast vertex). Such a block is called a “standard block”.

Affleck Oil Field, United Kingdom is an example of such a block.

2. The Needs of Transformation

Designing blocks is daily work in the resource exploration industry. Blocks serve many purposes, including but not limited to administration (relinquishing part of a block, dissolving several blocks into one), block bid campaigns, and so on.

Geologists usually draft blocks on a projected geological setting basemap. Unfortunately, the mapping software in use has no snap functionality. It can only export polygons (the drafted blocks) in projected XY format. It can, however, display a geodetic reference grid at 15-second intervals.

A new workflow was designed to handle the blocks:

1. Geologists design blocks whose vertices are visually, but not precisely, on the nodes of a 15-second geodetic reference grid.
2. All blocks reside in a layer, with valid identifiers.
3. Export the design layer of blocks in projected format as a text file.
4. Write a program to read the file and generate standard blocks.

3. Implementation

The implementation is a Python package called “transco” (TRANSformation of COordinates), written in 2022 and in use since then.

transco uses pyproj to perform the transformation. pyproj is a Python interface to PROJ. Shapely is also used. Shapely is a Python wrapper for GEOS.

3.1 class

transco uses three classes to represent the data: Layer, Block and Point. Points are properties of a Block instance, and blocks are properties of a Layer instance.

The Layer class is defined as follows:

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

class Layer:
    """
    represents a DoubleFox layer contains designed blocks.
    """

    def __init__(
        self,
        layer_file,
        coor_type="proj",
        crs="wgs84utm",
        location="location",
        compute_name=False,
    ):
        self.raw_lines = []
        self.bad_lines = []
        self._blocks = []
        self.cblocks = []
        self.coor_type = coor_type
        self.compute_name = compute_name
        self.name = layer_file.stem.replace("/", "_").replace("\\", "_")
        self.crs = get_crs(coor_type, crs, location)
        self.transformer = None
        self.transform_direction = "proj2geo"
        self.location = location
        self.transformer = get_transformer(self.crs, direction="proj2geo")
        with open(layer_file) as f:
            self.raw_lines = f.readlines()
        lines = copy.deepcopy(self.raw_lines)
        index_for_del = []
        for index, l in enumerate(lines):
            try:
                n, x, y = l.split(",")[:3]
                float(x)
                float(y)
            except:
                index_for_del.append(index)
                self.bad_lines.append(str(index + 1) + "," + l)
                continue
            finally:
                if not n:
                    raise Exception("Found point has no name")
        for i in sorted(index_for_del, reverse=True):
            del lines[i]
        groups = groupby(lines, key=lambda x: x.split(",")[0])
        self.raw_blocks = [
            [p.strip().split(",")[:3] for p in lines] for name, lines in groups
        ]
        for b in self.raw_blocks:
            block = Block(self, b)
            self._blocks.append(block)
        sorted_blocks = sorted(self.blocks, key=lambda x: x.name)
        block_groups = groupby(sorted_blocks, key=lambda x: x.name)
        raw_cblocks = [[b for b in cblock] for name, cblock in block_groups]
        for cb in raw_cblocks:
            _cb = CBlock(self, cb)
            self.cblocks.append(_cb)

    @property
    def blocks(self):
        if len(self._blocks) == len(self.raw_blocks):
            return self._blocks
        raise

    def __str__(self):
        return self.name + ":{} blocks".format(len(self.blocks))

    def __repr__(self):
        return self.__str__()

The coor_type, crs and location properties are used to obtain the corresponding CRS object from pyproj.

The Block class:

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

class Block:
    """A raw block resides in a Layer object."""

    def __init__(self, layer, block):
        self.layer = layer
        self.location = self.layer.location
        self.order = 0
        self.crs = self.layer.crs
        self.transformer = self.layer.transformer
        name_str = block[0][0]
        try:
            name, order = name_str.split("--")
            order = int(order)
            if name.strip():
                self.name = name.strip().replace("/", "").replace("\\", "")
            if order >= 1:
                self.order = order
        except:
            self.name = name_str.rstrip("--").strip().replace("/", "").replace("\\", "")
        if not self.name:
            raise Exception("Failed to create a name for this Block")
        self._points = []
        for p in block:
            point = Point(self, p)
            self._points.append(point)
        self.__regular_points()
        self.__unparallel_points()
        self.get_code_name()

    @property
    def points(self):
        return self._points

The block object gets CRS information from its Layer object.

And the Point class looks like:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

class Point:
    """represents a point resides in a Block object."""

    def __init__(self, block, point):
        self.block = block
        self.name = point[0].replace("/", "").replace("\\", "")
        self.x = float(point[1])
        self.y = float(point[2])
        self.lo_deg = None
        self.la_deg = None
        self.lo_dms = None
        self.la_dms = None
        self.offset_lo_dms = None
        self.offset_la_dms = None
        self.crs = self.block.crs
        self.transformer = self.block.transformer
        self._deg()
        self._dms()
        self._offset()
        self._offset_2_dms()

As before, a Point instance gets its CRS from its Block object.

3.2 Function

The basic idea is that the original XY coordinates are transformed into degrees and then snapped to the nearest geodetic reference grid node by the _offset method.

Several methods modify the blocks: rotating points if they are not clockwise, removing redundant points, checking for errors, computing block names by spatial relationship, and so on.

Here is the function for getting CRS.

def get_crs(coor_type, crs, location):
    try:
        epsg_code = defined_crs[coor_type][crs][location]
        return CRS.from_epsg(epsg_code)
    except:
        return None

Transform Degree to DMS:

1
2
3
4
5
6
7
8
9
10
11
12
13

def d2dms(x):
    """
    Transform a single float degree value to DMS format.
    return a tuple of degree, minute, second in float.
    >>> d2dms(113.5)
    (113.0, 30.0, 0.0)
    >>> d2dms(20.254166667)
    (20.0, 15.0, 15.000001199999247)
    """
    d_decimal, d = math.modf(x)
    m_decimal, m = math.modf(d_decimal * 60)
    s = m_decimal * 60
    return d, m, s

The last step is defining functions to call methods on the objects.

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

def load_file(
    file="df_design.txt", coor_type="proj", crs="wgs84utm", location="location"
):
    f = Path(file)
    return Layer(f, coor_type, crs, location)

def save_default(layer):
    layer.save_standardized_blocks()
    layer.save_blocks_area_info_wgs84utm()

def available_crs(defined_crs):
    all_crs = {}
    for key, value in defined_crs.items():
        if value:
            all_crs[key] = value
    return all_crs

def show_help(**kw):
    print(__desc__)
    print("Version {}".format(__version__))
    print("Author: ", __author__)
    print(__transco_help__)
    location = kw.get("location", None)
    if location:
        try:
            crs = get_crs("proj", "wgs84utm", location)
        except:
            print("Failed to get CRS.")
            raise
        if crs:
            print("Current CRS:")
            print(crs.name)
            print("Area of Use:")
            print(crs.area_of_use)

def main():
    parser = argparse.ArgumentParser(description=__desc__)
    parser.add_argument(
        "file", help="the input file, default is 'df_design.txt'.", nargs="?"
    )
    parser.add_argument(
        "location",
        help="the location for projection selection, default is 'location'.",
        nargs="?",
    )
    parser.add_argument(
        "-a", "--all", action="store_true", help="save all output files"
    )
    parser.add_argument("--version", action="store_true", help="show version")
    parser.add_argument("--allcrs", action="store_true", help="show all available CRS")
    args = parser.parse_args()
    try:
...

After that, I can call the main function from the CLI, call the functions and classes directly in an interactive Python environment, or use it in a Jupyter notebook.

4. Results

transco has been in use since 2022. A 50-block batch that previously took roughly an hour per block to prepare by hand in Excel now runs in about 10 minutes.

The larger gain is in the revision loop. Geologists revise their block designs repeatedly before a final version is agreed, and the manual cost was paid again on every revision. Re-running transco costs nothing, so a design can be iterated as many times as the work requires.

I tried to use Linux back in 2004, and started to use it as a daily operating system in 2016. Maybe because I wasn’t (though I tried to be one) an IT worker, I never met a person using Linux in person (using Redhat Workstation for Seismic Data Interpretation doesn’t count). I was asked by many people why using Linux?

Let me figure out an answer.

OpenWrt is an open source linux system tailored for routers. Using an offical OpenWrt VM system as router has obvious benefits: enjoy the most flexibility of OpenWrt and get rid of the outdated vendor firmware.

I have used my network scanner through Windows for a long time. It works well except I need to get a Windows laptop. Although the default storage path of the scanner resides in a samba server, it is still tedious to do daily scan. I managed to set it up on Linux. This blog is for my reference.