<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <a href="https
Step1:  Sample data
Step3:  A function to calculate three parameter transformation based on commmon points. The point coordinates are stored in dictionaries, the key is the point ID/name and each dictionary item stores a list of coordinates [x, y, z].
Step5:  Function to select two point from all points.
Step7:  Apply the transformation parameters to points.
Step8:  Iterating the transformation using two points in all combinations.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from math import atan2, sqrt, sin, cos, pi
import re
X, Y, Z, MX, MY, MZ = 0, 1, 2, 3, 4, 5   # indices in coord dictionary items
RO = 180 * 3600 / pi
Explanation: <a href="https://colab.research.google.com/github/OSGeoLabBp/tutorials/blob/master/english/data_processing/lessons/trans.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
Robust 2D transformation
Sometimes we have to check the stability or the horizontal reference system, which is established by the control points marked on the field.
The proposed method to check the movement of the control point:
Initially make observations and adjust the horizontal network as a free network with blunder elimination
Checking the stability of the control points, make observations and adjust the network again as a free network with blunder elimination
Calculate transformation parameters between the two free network using a robust method. If we suppose some control points moved between the two measurements, the moving points cannot be used to calculate the transformation parameters.
The robust method for the calculation of the transformation parameters can be the L1 norm, where we minimize the sum of the absolute value of corrections. It can be solved ba simplex method of linear programming or iterating the LSM solution using two points only for the parameter calculation. As the scale change between the two determinations of the points is not excepted 3 parameters ortogonal transformation (offset and rotation) is used.
Advanteges of the proposed method:
It is not neccessary to have the same points, there may be new and destroyed points
End of explanation
# initial coordinates
coo1 = {'K1': [ 0.0,     5.9427, 0.9950],
        'K2': [ 6.0242,  0.0,    1.3998],
        'K3': [ 9.7954,  5.3061, 1.8230],
        'K4': [17.9716,  5.2726, 1.8389],
        'K5': [31.6363,  5.5274, 1.0126],
        'K6': [33.2002,  7.0923, 1.1090],
        'K7': [35.9246, 14.5219, 1.3326],
        'K8': [40.6884, 21.0337, 1.4709],
        'K9': [32.501, 22.8658, 1.6797]
       }
# coordinatess from the second adjustment
coo2 = {'K1': [ 0.0002,  5.9422, 0.9948],
        'K2': [ 6.0252, -0.0006, 1.3997],
        'K3': [ 9.7959,  5.3061, 1.8230],
        'K4': [17.9716,  5.2729, 1.8389],
        'K5': [31.6366,  5.5280, 1.0129],
        'K6': [33.1994,  7.0916, 1.1091],
        'K7': [35.9235, 14.5207, 1.3327],
        'K8': [40.6888, 21.0319, 1.4711],
        'K9': [32.2494, 22.8644, 1.6799]
       }
Explanation: Sample data
End of explanation
def tr3(src, dst, x0):
     Three parameter orthogonal transformation
        :param src: dictionary of source points and coordinates
        :param dst: dictionary of target points and coordinates
        :param x0: preliminary transformation parameter values
        :returns: x_offset y_offset rotation
    
    # find common points
    s = set(src.keys())
    d = set(dst.keys())
    common = s.intersection(d)
    n = len(common)
    A = np.zeros((2*n, 3))
    l = np.zeros(2*n)
    i = 0
    # set up equations
    for key in common:
        A[i] = np.array([1.0, 0.0, -src[key][X] * sin(x0[2]) -
                         src[key][Y] * cos(x0[2])])
        l[i] = dst[key][X] - (x0[0] + src[key][X] * cos(x0[2]) -
                              src[key][Y] * sin(x0[2]))
        i += 1
        A[i] = np.array([0.0, 1.0, src[key][X] * cos(x0[2]) -
                         src[key][Y] * sin(x0[2])])
        l[i] = dst[key][1] - (x0[1] + src[key][X] * sin(x0[2]) +
                              src[key][Y] * cos(x0[2]))
        i += 1
    # solve equation
    ATA = np.dot(A.transpose(), A)
    ATl = np.dot(A.transpose(), l)
    param = np.linalg.solve(ATA, ATl)   # x0, y0, rotation
    v = np.dot(A, param+x0) - l         # corrections
    return param + x0, v
Explanation: A function to calculate three parameter transformation based on commmon points. The point coordinates are stored in dictionaries, the key is the point ID/name and each dictionary item stores a list of coordinates [x, y, z].
End of explanation
def sel(coo, keys):
     select points from coordinate list based on point IDs or regexp
        :param coo: dictionary with coordinates
        :param keys: dictionary keys/point IDS to select or a regexp for point ids
    
    if isinstance(keys, str):
        r = re.compile(keys)
        w = list(filter(r.search, coo.keys()))
    else:
        w = keys
    return {k : coo[k] for k in w if k in coo}
Explanation: Function to select two point from all points.
End of explanation
def coo_tr(coo, param):
     transform coordinates in coo using transformation parameters
        :param coo: dictionary of coordinates to transform
        :param param: transformation parameters x0, y0, alfa, scale
    
    if len(param) == 4:
        x0, y0, alpha, scale = param
    else:
        x0, y0, alpha = param
        scale = 1.0
    return {k: [x0 + coo[k][X] * scale * cos(alpha) - coo[k][Y] * scale * sin(alpha),
                y0 + coo[k][X] * scale * sin(alpha) + coo[k][Y] * scale * cos(alpha),
                coo[k][Z]] for k in coo}
Explanation: Apply the transformation parameters to points.
End of explanation
key_list = list(coo1.keys())
n_key = len(key_list)
min_v = 1e38
print('P1   P2      X0       Y0     Alpha"   sum(|v|)')
print('----------------------------------------------')
for i in range(n_key):
    k1 = key_list[i]
    for j in range(i+1, n_key):
        k2 = key_list[j]
        p, v = tr3(sel(coo1, [k1, k2]), sel(coo2, [k1, k2]), [0.0, 0.0, 0.0])
        coo1_tr = coo_tr(coo1, p)
        sum_v = 0
        # calculate sum of absolute value of corrections
        for k in coo1:
            sum_v += abs(coo1_tr[k][X] - coo2[k][X]) + \
                     abs(coo1_tr[k][Y] - coo2[k][Y])
        if sum_v < min_v:
            opt = [k1, k2, p, sum_v]
            min_v = sum_v
        print(f'{k1:4s} {k2:4s} {p[0]:8.3f} {p[1]:8.3f} {p[2] * RO:6.1f} {sum_v:8.3f}')
print('optimal:')
print(f'{opt[0]:4s} {opt[1]:4s} {opt[2][0]:8.3f} {opt[2][1]:8.3f} {opt[2][2] * RO:6.1f} {opt[3]:8.3f}')
Explanation: Iterating the transformation using two points in all combinations.
End of explanation
<END_TASK>