Massimo di Stefano (aka Epifanio) is a guest blogger on Spatialguru.com

Feb 14 14:53

Datum transformation in Python

By: epifanio

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
»
Feb 10 00:50

Epsg code in python

By: epifanio

Python script to store epsg code, params and title in a dictionary and a simple search tool.

It's just a try ... episg.py

How to run in python : 

    >>from episg import *
 
    >>test()  # print out some examples
 
	#  rep3('path to epsg','order',str('order'),'a')
	#  order (code,param,title)
	# 'c' epsg code 
	# 't' title
	# 'p' parameters
	# 'a' all
	# input : epsg code
	# output param, code, title (all)
 
	>>>output = rep3('/Users/sasha/Desktop/epsg','code',str('4326'),'a')
	>>>print output
 
['proj=longlat', 'ellps=WGS84', 'datum=WGS84', 'no_defs'] 4326 WGS 84
How to run in bash : 
    $python episg.py /Users/sasha/Desktop/epsg code 3004 t
    $Monte Mario - Italy zone 2
The epsg file comes from the proj4-source code.