Forum

Author Topic: Translating WGS84 coordinates to local coordinates  (Read 404 times)

colin.spexigeo

  • Newbie
  • *
  • Posts: 5
    • View Profile
Translating WGS84 coordinates to local coordinates
« on: October 05, 2018, 09:15:14 PM »
Hi there,

I'm trying to convert a polygon with WGS84 coordinates to local space so I can fit Agisoft's region box around a predefined geofence.

I have tried the following, but It gives me an "Unsupported datum transformation" exception.

Code: [Select]
chunk.crs
Out[372]: 2018-10-05 11:11:56 <CoordinateSystem 'WGS 84 (EPSG::4326)'>

chunk.shapes.crs
Out[373]: 2018-10-05 11:12:10 <CoordinateSystem 'WGS 84 (EPSG::4326)'>

chunk.shapes[0].vertices
Out[374]: 2018-10-05 11:12:17 Shape.Vertices([Vector([-123.06583, 49.02038, 0.0]), Vector([-123.065345, 49.020376, 0.0]), Vector([-123.06535, 49.020779, 0.0]), Vector([-123.065825, 49.020783, 0.0]), Vector([-123.06583, 49.02038, 0.0])])

localCrs = PhotoScan.CoordinateSystem('LOCAL')

localCrs
Out[376]: 2018-10-05 11:13:43 <CoordinateSystem 'Local Coordinates (m)'>

PhotoScan.CoordinateSystem.transform(chunk.shapes[0].vertices[0], chunk.shapes.crs, localCrs)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-377-d48b461d2fab> in <module>()
----> 1 PhotoScan.CoordinateSystem.transform(chunk.shapes[0].vertices[0], chunk.shapes.crs, localCrs)

RuntimeError: Unsupported datum transformation

Any help is appreciated! Thank you!

Alexey Pasumansky

  • Agisoft Technical Support
  • Hero Member
  • *****
  • Posts: 9917
    • View Profile
Re: Translating WGS84 coordinates to local coordinates
« Reply #1 on: October 05, 2018, 10:13:39 PM »
Hello colin.spexigeo,

If you need to transform the coordinates from the shapes coordinate system to the internal coordinate system of the chunk (in which chunk.region.center, for example, is defined), then I can suggest the following code:

Code: [Select]
chunk = PhotoScan.app.document.chunk
T = chunk.transform.matrix
crs = chunk.shapes.crs

s = chunk.shapes[0]
coords = [T.inv().mulp(crs.unproject(v)) for v in s.vertices]
Best regards,
Alexey Pasumansky,
AgiSoft LLC

colin.spexigeo

  • Newbie
  • *
  • Posts: 5
    • View Profile
Re: Translating WGS84 coordinates to local coordinates
« Reply #2 on: October 09, 2018, 08:45:27 PM »
Thanks Alexey, works great!

colin.spexigeo

  • Newbie
  • *
  • Posts: 5
    • View Profile
Re: Translating WGS84 coordinates to local coordinates
« Reply #3 on: October 10, 2018, 11:53:48 PM »
For posterity and for those interested, I finally put all the pieces together on the thing I was working on. This script does the following:

1. For a given polygon shape, create a convex hull.
2. From the convex hull, determine the minimum-area oriented bounding box (obb)
3. Set the chunk region to be this oriented bounding box (conform_region_to_polygon).

This was all needed because we want to use known GEOJSON shapes of property lines to shrink the processing region automatically.

NOTE: the script assumes that:
1. You've already run "Align Photos..." on a set of georeferenced photos
2. The chunk coordinate system is WGS 84.
3. The chunk has a single polygon shape with at least 3 vertices.

Code: [Select]
__author__ = 'Colin Basnett (https://github.com/cmbasnett)'

import math
from functools import cmp_to_key
from PhotoScan import *


# Three points are a counter-clockwise turn if ccw > 0, clockwise if
# ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
# gives twice the signed  area of the triangle formed by p1, p2 and p3.
def ccw(p1, p2, p3):
    return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x)


def angle(p1, p2):
    x = (p2.x - p1.x)
    y = (p2.y - p1.y)
    a = math.atan2(y, x)
    return a


def points_cmp(a, b):
    if a.y == b.y:
        return a.x - b.x
    else:
        return a.y - b.y


# Removes duplicates from a sorted list (duplicate items are guaranteed to be adjacent)
def remove_duplicates(points):
    old = points[-1]
    for i in range(len(points) - 2, 0, -1):
        if points[i] == old:
            points.pop(i)
        else:
            old = points[i]


# Graham Scan to get a convex hull from a set of points
def get_convex_hull(points):
    # Sort by y, then x
    points = sorted(points, key=cmp_to_key(points_cmp))
    # Remove duplicates from sorted list
    remove_duplicates(points)
    # Now sort by polar angle from the first point
    points[1:] = sorted(points[1:], key=lambda p: angle(points[0], p))
    # Remove adjacent collinear line segments.
    for i in range(len(points) - 2, 0, -1):
        a = angle(points[i - 1], points[i])
        b = angle(points[i], points[i + 1])
        if a == b:
            points.pop(i)
    stack = points[0:3]
    for i in range(3, len(points)):
        while ccw(stack[-2], stack[-1], points[i]) <= 0:
            stack.pop()
        stack.append(points[i])
    return stack


def distance(p1, p2):
    return math.hypot(p1.x - p2.x, p1.y - p2.y)


# Numerically stable version of Heron's formula
def area_of_triangle(points):
    assert len(points) == 3
    a, b, c = points
    a, b, c = sorted([distance(a, b), distance(b, c), distance(c, a)], reverse=True)
    return 0.25 * math.sqrt((a + (b + c)) * (c - (a - b)) * (c + (a -b)) * (a + (b - c)))


# Area of polygon as a sum of triangle areas
def area_of_convex_hull(convex_hull):
    poly_area = 0
    for i in range(1, len(convex_hull) - 1, 1):
        points = [convex_hull[0], convex_hull[i], convex_hull[i + 1]]
        poly_area += area_of_triangle(points)
    return poly_area


def rotate_points(points, origin, angle):
    new_points = list(map(lambda x: x - origin, points))
    for p in new_points:
        s = math.sin(angle)
        c = math.cos(angle)
        px = (p.x * c) - (p.y * s)
        py = (p.x * s) + (p.y * c)
        p.x = px
        p.y = py
    return list(map(lambda x: x + origin, new_points))


def fit_obb(points):
    min_area = None
    obb = None
    for i in range(len(points) - 1):
        theta = -angle(points[i], points[i + 1])
        origin = points[i]
        rotated_polygon = rotate_points(points, origin, theta)
        aabb = aabb_from_points(rotated_polygon)
        area = aabb.area()
        if min_area is None or area < min_area:
            min_area = area
            size = aabb.size()
            rot = -theta
            obb_vertices = rotate_points(aabb.vertices(), origin, rot)
            center = (obb_vertices[0] + obb_vertices[2]) * 0.5
            obb = OBB(center, rot, size)
    return obb


class OBB(object):
    def __init__(self, center, rot, size):
        self.center = center
        self.rot = rot
        self.size = size

    def vertices(self):
        points = [
            Vector([-self.size.x / 2, -self.size.y / 2, 0]),
            Vector([ self.size.x / 2, -self.size.y / 2, 0]),
            Vector([ self.size.x / 2, self.size.y / 2, 0]),
            Vector([-self.size.x / 2, self.size.y / 2, 0]),
        ]
        points = rotate_points(points, Vector([0, 0, 0]), self.rot)
        return list(map(lambda x: x + self.center, points))


def aabb_from_points(points):
    min_x = min(points, key=lambda p: p.x).x
    min_y = min(points, key=lambda p: p.y).y
    max_x = max(points, key=lambda p: p.x).x
    max_y = max(points, key=lambda p: p.y).y
    return AABB(Vector([min_x, min_y, 0]), Vector([max_x, max_y, 0]))


class AABB(object):
    def __init__(self, min_, max_):
        self.min = min_
        self.max = max_

    def center(self):
        return (self.max + self.min) * 0.5

    def area(self):
        a = self.max - self.min
        return a.x * a.y

    def vertices(self):
        return [
            self.min,
            Vector([self.max.x, self.min.y, 0]),
            self.max,
            Vector([self.min.x, self.max.y, 0])
        ]

    def size(self):
        return self.max - self.min


def set_rot_of_matrix(transform, rot):
    transform[0, 0] = rot[0, 0]
    transform[0, 1] = rot[0, 1]
    transform[0, 2] = rot[0, 2]
    transform[1, 0] = rot[1, 0]
    transform[1, 1] = rot[1, 1]
    transform[1, 2] = rot[1, 2]
    transform[2, 0] = rot[2, 0]
    transform[2, 1] = rot[2, 1]
    transform[2, 2] = rot[2, 2]


def region_transform(chunk):
    transform = Matrix.Translation(chunk.region.center)
    set_rot_of_matrix(transform, chunk.region.rot)
    return transform


def local_to_region(chunk, points):
    transform = region_transform(chunk)
    points = [transform.inv().mulp(x) for x in points]
    return points


def convert_to_region(chunk, points):
    T = chunk.transform.matrix
    crs = chunk.shapes.crs
    points = [T.inv().mulp(crs.unproject(x)) for x in points]
    return local_to_region(chunk, points)


# Shrinks the region to a minimum bounding box around a polygon.
def conform_region_to_polygon(chunk, shape):
    assert chunk is not None
    assert shape is not None
    assert shape.type is Shape.Polygon
    assert len(shape.vertices) >= 3
    chunk.resetRegion()
    # Save the old region rotation as we will use it in the final region translation.
    old_region_rot = chunk.region.rot
    # Get the WSG vertices in region space
    vertices = convert_to_region(chunk, chunk.shapes[0].vertices)
    # Generate a convex hull around the shape vertices
    convex_hull = get_convex_hull(vertices)
    # Fit a minimal bounding box around the convex hull
    obb = fit_obb(convex_hull)
    # Maintain the chunk region height
    obb.size.z = chunk.region.size.z
    # Resize the chunk to the size of the OBB
    chunk.region.size = obb.size
    # Rotate the region by the rotation of the OBB
    chunk.region.rot = chunk.region.rot * Utils.ypr2mat(Vector([math.degrees(obb.rot), 0, 0])).inv()
    # Translate the region center to the OBB's center
    transform = Matrix.Diag([1, 1, 1, 1])
    set_rot_of_matrix(transform, old_region_rot)
    chunk.region.center += transform.mulp(obb.center)


# Do the deed!
chunk = PhotoScan.app.document.chunk
conform_region_to_polygon(chunk, chunk.shapes[0])