Recent blog posts
Guest Blogs
Syndicate
Datum transformation in Python
I’ll describe here a procedure for the evaluation of the tranforming parameters necessary to effect coordinates change among reference frames materialized on different Datum.
In order to have the changing of Datum it is used a geometric transformation of correspond (conforme) type through the resolution of a last squares system based on a homologous points series, in order toobtain the direct changing from plane coordinates in one of the two Datum to the corresponding in the other Datum.
The procedure realyzed in Python environment, provides the possibilty to apply both a conforme transformation (plane rotation with translation and isotropic scale variation) and of affine transformation (plane rotation with translation and anisotropic scale variation).
The code has been implemented in a Python module sasha.py, in order to operate geometric transformation on text files, trhough the use of special software libraries (numpy , scipy ) .
An example on how to use scipy-numpy to make a datum transformation :>>from numpy import *
>>from scipy import *
>>L = io.readarray(str('locale'))
>>G = io.readarray(str('globale'))
>>L
array([[ 496.4 , -14853.55],
[ -405.58, -15197.24],
[ 602.12, -15273.71],
[ 643.14, -15141.68]])
>>G
array([[ 2561715.65, 4444638.81],
[ 2560815.52, 4444289.81],
[ 2561823.58, 4444219.83],
[ 2561863.75, 4444351.03]])
>>A = zeros((2*L.shape[0],4),float)
>>A
array([[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.]])
>>A[ ::2, 0] = 1.0
>>A[1::2, 1] = 1.0
>>A[ ::2, 2] = L[:,0]
>>A[1::2, 2] = L[:,1]
>>A[ ::2, 3] = L[:,1]
>>A[1::2, 3] = -L[:,0]
>>A
array([[ 1.00000000e+00, 0.00000000e+00, 4.96400000e+02,
-1.48535500e+04],
[ 0.00000000e+00, 1.00000000e+00, -1.48535500e+04,
-4.96400000e+02],
[ 1.00000000e+00, 0.00000000e+00, -4.05580000e+02,
-1.51972400e+04],
[ 0.00000000e+00, 1.00000000e+00, -1.51972400e+04,
4.05580000e+02],
[ 1.00000000e+00, 0.00000000e+00, 6.02120000e+02,
-1.52737100e+04],
[ 0.00000000e+00, 1.00000000e+00, -1.52737100e+04,
-6.02120000e+02],
[ 1.00000000e+00, 0.00000000e+00, 6.43140000e+02,
-1.51416800e+04],
[ 0.00000000e+00, 1.00000000e+00, -1.51416800e+04,
-6.43140000e+02]])
>>Y = zeros((2*G.shape[0],1),float)
array([[ 0.],
[ 0.],
[ 0.],
[ 0.],
[ 0.],
[ 0.],
[ 0.],
[ 0.]])
>>Y[ ::2, 0] = G[:,0]
>>Y[1::2, 0] = G[:,1]
>>Y
array([[ 2561715.65],
[ 4444638.81],
[ 2560815.52],
[ 4444289.81],
[ 2561823.58],
[ 4444219.83],
[ 2561863.75],
[ 4444351.03]])
>>N = dot(A.T.conj(), A)
>>T = dot(A.T.conj(), Y)
>>C = dot(linalg.inv(N), T)
>>C
array([[ 2.56113298e+06],
[ 4.45948730e+06],
[ 9.99856124e-01],
[ -5.80009571e-03]])
>>E0 = C[0]
>>N0 = C[1]
>>a = C[2]
>>b = C[3]
>>x = E0 + L[0,0] * a + L[0,1] * b
>>y = N0 + L[0,1] * a - L[0,0] * b
>>print x,y
[ 2561715.45624184] [ 4444638.76898143]Examples on how to use sasha.py
>> #in python
>> from sasha import conforme
>> print conforme('GCP1.txt','GCP2.txt','x','out.txt')
Parametri di trasformazione :
Tx : -2587285.48816
Ty : -4445117.1582
Rotazione rigida : [ 0.33236506]
Fattore di scala: [ 1.]
Varianza dei parametri geometrici :
Sigma quadro rotazione rigida : 0.0244664334334
Sigma quadro varazione di scala : 0.000854147546627
Scarto quadratico medio : 0.343476930718
none
Post new comment