Fuente de obtención de datos

Para la demostración de la obtención de los ángulos de rotación en el proceso de orientación interna, se seleccionó una imagen de fotografía aérea desde el repositorio de Colombia en mapas que se ubica en la siguiente url: clic aquí , en seguida, se descargó la propia imagen,obtenido del sensor Wild RC 30 y el certificado de calibración, para obtener los valores de las coordenadas de las marcas fiduciales en la propia imagen y en el certificado de calibración como se observa en la siguiente figura:

Datos de coordenadas de las marcas fiduciales de la imagen y del certificado de calibración
Datos de coordenadas de las marcas fiduciales de la imagen y del certificado de calibración

Consideraciones teóricas

Teniendo en cuenta que en orientación interior, ocurre el proceso de transformación bidimensional de tipo afín, donde participa el ángulo de giro, el ángulo de compensación de perpendicularidad, el factor de escala en el eje X y el actor de escala en el aje Y, los valores de traslación en el eje “X” y “Y”, para lograr la orientación interior una las opciones más recomendadas es aplicando la siguiente expresión:

\[\begin{equation} \begin{bmatrix} x\\ y \end{bmatrix}= \begin{bmatrix} cos\alpha &-sen\alpha\\ sen\alpha & cos\alpha \end{bmatrix}. \begin{bmatrix} 1 &-sen\beta\\ 0 & cos\beta \end{bmatrix}. \begin{bmatrix} \lambda_{x}x^{'}\\ \lambda_{y}y^{'} \end{bmatrix}+ \begin{bmatrix} T_{x}\\ T_{y} \end{bmatrix} \end{equation}\]

En ese orden de ideas, previamente es posible simplificar la expresión anterior, quedando la ecuación como sigue:

\[\begin{equation} \begin{bmatrix} X\\ Y \end{bmatrix}= \begin{bmatrix} a&b\\ c&d \end{bmatrix}. \begin{bmatrix} x\\ y \end{bmatrix}+ \begin{bmatrix} T_{x}\\ T_{y} \end{bmatrix} \end{equation}\]

Donde la intención es conocer los valores de: “a”, “b” “c” “d”, \(T_{x}\) y \(T_{y}\), a partir ecuaciones lineales de seis incógnitas, para dar solución de este inconveniente, nos aproximamos a resolver utilizando el criterios de mínimos cuadrados; para mayor detalle, revisar el manuscrito de apuntes de fotogrametría II se se puede acceder en la siguiente url: acceder aquí y revisar la pagina 25.

Sistema de ecuación matricial, fuente de imagen (Perez, 2001)
Sistema de ecuación matricial, fuente de imagen (Perez, 2001)

Al resolver la ecuación se obtiene los valores de “a”, “b” “c” “d”, \(T_{x}\) y \(T_{y}\), los valores de los residuales, el valor de la precisión de la transformación, este ultimo aplicando la siguiente expresión:

\[\begin{equation} \sigma_{xy} = \sqrt{\frac{\sum_{i=1}^{n} (V^{2}_{xi}+V^{2}_{yj})}{l-m}} \end{equation}\]

donde:

l = “es el número de relaciones de observación (número de ecuaciones) que intervienen en el proceso de cálculo de los parámetros, es decir, l = 2n”, (Perez, 2001).

m= “es el número de relaciones de observación mínimas que se necesitan para resolver los parámetros, en esta transformación m=4”, (Perez, 2001).

Así mismo, de los resultados, se posible conocer los valores de factor de escala, para los ejes X y Y, el valor del ángulo de rotación (\(\alpha\)), el valor del ángulo de compensación de perpendicularidad (\(\beta\)), utilizando las siguientes relaciones:

Cálculo de ángulo de rotación.

\[\begin{equation} \alpha = atan(\frac{c}{a}) \end{equation}\]

Cálculo de ángulo de compensación de perpendicularidad.

\[\begin{equation} \beta = \alpha +atan(\frac{b}{d}) \end{equation}\]

Cálculo de valor de escala en el eje X.

\[\begin{equation} \lambda_{x}= \frac{a}{cos\alpha} \end{equation}\]

Cálculo de valor de escala en el eje Y.

\[\begin{equation} \lambda_{y}= \frac{d}{cos(\alpha-\beta)} \end{equation}\]

Ejercicio Aplicativo

Importando los datos.

#dato <-  read_excel("D:/.../DatoOIm.xlsx", sheet = "A")

A <- data.matrix(dato)
A
##               1         2         3         4 5 6
##  [1,] 0.2203891 0.0073742 0.0000000 0.0000000 1 0
##  [2,] 0.0000000 0.0000000 0.2203891 0.0073742 0 1
##  [3,] 0.0083581 0.0079342 0.0000000 0.0000000 1 0
##  [4,] 0.0000000 0.0000000 0.0083581 0.0079342 0 1
##  [5,] 0.0088878 0.2199374 0.0000000 0.0000000 1 0
##  [6,] 0.0000000 0.0000000 0.0088878 0.2199374 0 1
##  [7,] 0.2209024 0.2193982 0.0000000 0.0000000 1 0
##  [8,] 0.0000000 0.0000000 0.2209024 0.2193982 0 1
##  [9,] 0.1143544 0.0016628 0.0000000 0.0000000 1 0
## [10,] 0.0000000 0.0000000 0.1143544 0.0016628 0 1
## [11,] 0.0026287 0.1139587 0.0000000 0.0000000 1 0
## [12,] 0.0000000 0.0000000 0.0026287 0.1139587 0 1
## [13,] 0.1149132 0.2256663 0.0000000 0.0000000 1 0
## [14,] 0.0000000 0.0000000 0.1149132 0.2256663 0 1
## [15,] 0.2266391 0.1133667 0.0000000 0.0000000 1 0
## [16,] 0.0000000 0.0000000 0.2266391 0.1133667 0 1
# dato1 <-  read_excel("D:/.../DatoOIm.xlsx", sheet = "L")

L <-data.matrix(dato1)
L
##              1
##  [1,]  106.000
##  [2,] -105.995
##  [3,] -105.999
##  [4,] -105.999
##  [5,] -105.998
##  [6,]  105.993
##  [7,]  106.000
##  [8,]  106.000
##  [9,]   -0.006
## [10,] -111.999
## [11,] -111.999
## [12,]   -0.002
## [13,]    0.004
## [14,]  111.998
## [15,]  112.005
## [16,]   -0.001
Aplicando las ecuaciones
# Hallando los valores de las variables no conocidas.
X<- (solve(t(A)%*%A,tol = 1e-200))%*%(t(A)%*%L)
X
##             1
## 1  999.909802
## 2   -2.454506
## 3    2.628331
## 4  999.922261
## 5 -114.343900
## 6 -113.955398
X[4,1]
## [1] 999.9223
# Cálculo del ángulo  de giro α 

α <- atan(X[3,1]/X[1,1])*(180/pi)
α
## [1] 0.1506055
# Cálculo del ángulo  de falta de perpendicularidad  ß 

ß <- atan(X[2,1]/X[4,1])*(180/pi)+α
ß
## [1] 0.009962068
#Cálculo del factor de escala en el eje X 
Ex <- X[1,1]/cos((α*pi)/180)
Ex
## [1] 999.9133
#Cálculo del factor de escala en el eje X 
Ey <- X[4,1]/cos(((α-ß)*pi)/180)
Ey
## [1] 999.9253
V = A%*%X-L
V
##                   1
##  [1,]  7.220907e-03
##  [2,] -7.515589e-03
##  [3,] -7.028882e-03
##  [4,] -8.468644e-04
##  [5,]  1.260279e-03
##  [6,] -4.735532e-03
##  [7,]  6.049313e-05
##  [8,]  6.351022e-03
##  [9,]  2.103669e-03
## [10,]  6.834074e-03
## [11,]  3.850158e-03
## [12,]  3.352146e-03
## [13,]  1.035401e-03
## [14,] -2.610999e-03
## [15,] -8.502024e-03
## [16,] -8.282573e-04
# Cálculo de la precisión de la transformación el valor de 8 es número de puntos de control,
# el valor de 6 es número mínimo observaciones a resolver. 

s <- (sum(V^2)/(2*8-6))^2
s
## [1] 1.452761e-09
### Construyendo el modelo

oi<- function(x=0.21469 ,y=0.19731,α1=α*pi/180, ß1=ß*pi/180,tx=X[5,1], ty=X[6,1],E1=Ex,E2=Ey){
  es <- matrix(c(E1,0,0,E2),ncol =2,byrow = TRUE)
  rt <- matrix(c(cos(α1),-sin(α1),sin(α1),cos(α1)),ncol = 2,byrow = TRUE)
  cp <- matrix(c(1,-sin(ß1),0,cos(ß1)),nrow = 2,byrow = TRUE)
  ta <- matrix(c(tx,ty),ncol =1,byrow = TRUE)
  d <-  matrix(c(x,y),ncol =1,byrow = TRUE)
  p <-  es%*%rt%*%cp%*%d+ta
  return(p)
}
oi()
##          [,1]
## [1,] 99.77384
## [2,] 83.90337
plot(oi()[1,1],oi()[2,1])

## Warning: package 'sp' was built under R version 4.3.3
##  [1] "C-2748-0024.mgrd"                       
##  [2] "C-2748-0024.prj"                        
##  [3] "C-2748-0024.sdat"                       
##  [4] "C-2748-0024.sgrd"                       
##  [5] "Codigo  transformacion de datos.txt"    
##  [6] "Codigo modelo transformacion 3D.txt"    
##  [7] "Codigo Transformacion bidimencional.txt"
##  [8] "data.csv"                               
##  [9] "data.txt"                               
## [10] "raster.mgrd"                            
## [11] "raster.prj"                             
## [12] "raster.sdat"                            
## [13] "raster.sgrd"
r
## class      : RasterLayer 
## dimensions : 2584, 4043, 10447112  (nrow, ncol, ncell)
## resolution : 1.5003e-05, 1.5003e-05  (x, y)
## extent     : 0.1234147, 0.1840718, 0.04353871, 0.08230646  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C-2748-0024.sdat 
## names      : C.2748.0024
plot(r,col = grey.colors(255))

# Extrayendo las coordenadas de los pixeles del dato raster.
a <- xFromCell(r,1:10447112)

b<- yFromCell(r,1:10447112)

Aplicando la ecuación de la transformación afín (primera ecuación), y remplazando los valores calculados de : \(\alpha\), \(\beta\), \(\lambda_{x}\) = Ex, \(\lambda_{y}\) = Ey, \(T_{x}\) y \(T_{y}\).

## class      : RasterLayer 
## dimensions : 2584, 4043, 10447112  (nrow, ncol, ncell)
## resolution : 0.0150248, 0.0150577  (x, y)
## extent     : 8.836533, 69.58179, -70.08822, -31.17913  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 24, 249  (min, max)
roi <-m()
roi
## class      : RasterLayer 
## dimensions : 2584, 4043, 10447112  (nrow, ncol, ncell)
## resolution : 0.0150248, 0.0150577  (x, y)
## extent     : 8.836533, 69.58179, -70.08822, -31.17913  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 24, 249  (min, max)
# Evidenciando el proceso de orientación interior 
par(mfrow = c(1,1))
plot(roi,col = grey.colors(255))

#Comparando los productos antes del proceso de transformación y después de la transformación. 
par(mfrow = c(2,2))
plot(r,col = grey.colors(255))
plot(roi,col = grey.colors(255))