Script para georreferenciar archivos vectoriales, con puntos comunes

vectorwrap.py

#!/usr/bin/env python3

import numpy as np
from pandas import read_table
from io import StringIO,BytesIO
import logging
import subprocess as sp
import click

log = logging.getLogger()

def readGCPs(gcps):
    """
    """
    log.debug(gcps)
    try:
        table = read_table(gcps,sep=None,header=None,dtype=float)
        return table
    except ValueError:
        try:
            table = read_table(gcps,sep=None,dtype=float)
            return table
        except ValueError:
            log.debug('not a file or file-like')

    try:
        sio = StringIO(gcps)
    except TypeError:
        log.debug('not a string')
        sio = StringIO('\n'.join(gcps))
        #Any Exception here will pass uncatched

    try:
        table = read_table(sio,sep=None,header=None,dtype=float)
        return table
    except ValueError:
        table = read_table(sio,sep=None,dtype=float)
        return table
        #Any exception Here will pass uncatched

def kabsch_umeyama(A, B):
    """
    Fuente:
    https://zpl.fi/aligning-point-patterns-with-kabsch-umeyama-algorithm/
    """
    assert A.shape == B.shape
    n, m = A.shape

    EA = np.mean(A, axis=0)
    EB = np.mean(B, axis=0)
    VarA = np.mean(np.linalg.norm(A - EA, axis=1) ** 2)

    H = ((A - EA).T @ (B - EB)) / n
    U, D, VT = np.linalg.svd(H)
    d = np.sign(np.linalg.det(U) * np.linalg.det(VT))
    S = np.diag([1] * (m - 1) + [d])

    R = U @ S @ VT
    c = VarA / np.trace(np.diag(D) @ S)
    t = EA - c * R @ EB

    return R, c, t

def getTransform(gcp_table,method='Procustes'):
    """
    """

    in_points = gcp_table.values[:,:2]
    #in_points = np.stack([*in_points.T,np.zeros(len(in_points))],axis=1)
    out_points = gcp_table.values[:,2:]
    #out_points = np.stack([*out_points.T,np.zeros(len(out_points))],axis=1)

    print(in_points)
    print(out_points)

    R,c,t = kabsch_umeyama(out_points,in_points)

    print(R)
    print(c)
    print(t)

    ((s11,s12),(s21,s22)) = c * R
    (xoff,yoff) = t

    affine_step = ['+proj=affine',f'+s11={s11}',f'+s12={s12}'
                                 ,f'+s21={s21}',f'+s22={s22}'
                                 ,f'+xoff={xoff}',f'+yoff={yoff}']

    return affine_step

def transformGCPs(gcp_table,transform):
    in_points = gcp_table.values[:,:2]
    command = ['cct','-z0','-t0',*transform.split()]

    bytes_in = '\n'.join(map(lambda x: '{} {}'.format(*x),in_points)).encode()+b'\n'
    new_in_points = sp.check_output(command,input=bytes_in)
    bio = BytesIO(new_in_point)
    table = read_table(bio,sep=None,header=None,dtype=float)
    new_gcp_table = gcp_table.copy()
    new_gcp_table.loc[:,:2] = table.loc[:,:2]

    return new_gcp_table

@click.command()
@click.argument('vector_in',type=str)
@click.argument('vector_out',type=str)
@click.option('--gcps',type=str,multiple=True)
@click.option('--gcp_transform',type=str)
@click.option('--ogroptions',type=str,multiple=True)
def vectorWrap(vector_in, vector_out, gcps, gcp_transform=None, ogroptions=[]):
    gcps = list(gcps)

    if len(gcps) == 1:
        gcps = gcps[0]

    gcp_table = readGCPs(gcps)

    if not gcp_transform is None:
        gcp_table = transformGCPs(gcp_table,gcp_transform)

    mapping = getTransform(gcp_table)

    command = ['ogr2ogr','-ct',' '.join(mapping),*ogroptions,vector_out,vector_in]

    p = sp.run(command)

    return p.returncode

if __name__ == '__main__':
    log.setLevel(logging.DEBUG)
    exit(vectorWrap())