En este post doy un ejemplo de cómo usar mi primer paquete de R, llamado ‘vec2dtransf‘, para realizar transformaciones afínes y de semejanza sobre datos vectoriales.
For English click here.
Casos de uso
Las transformaciones afínes y de semejanza son útiles cuando se integran datos espaciales de distintas fuentes. Suele ocurrir que los vectores de un conjunto de datos (llamémoslo ‘A’) no coinciden con un conjunto de datos base (‘B’), el cual puede ser raster o vector. En dicho escenario uno buscaría reposicionar ‘A’ tomando como referencia ‘B’.
Hay varios casos en los que esta situación ocurre, por ejemplo:
- ‘A’ no tiene información de su proyección y no hay ninguna pista que la indique.
- ‘A’ ya fue proyectado pero simplemente no coincide con ‘B’ (piensen en datos con datum local cuando se migran a WGS84).
- ‘A’ tiene un sistema de coordenadas arbitrario o falso y va a ser integrado con datos de carácter nacional. Este puede ser el caso de proyectos arqueológicos o de ingeniería.
- ‘A’ fue digitalizado a partir de una imagen distorsionada, como un mapa antiguo escaneado.
- ‘A’ fue digitalizado a partir de una imagen de satélite mal georreferenciada (recientemente ocurrió en OpenStreetMap).
El paquete ‘vec2dtransf’
El paquete de R ‘vec2dtransf’ provee clases para definir y realizar transformaciones afínes y de semejanza sobre datos vectoriales, concretamente en objetos ‘sp’ (objetos creados por el paquete de R ‘sp‘). Las transformaciones se pueden definir a partir de puntos de control o directamente desde parámetros de transformación. Si se proveen más puntos de control que los requeridos, se ejecutan mínimos cuadrados, permitiendo obtener residuales y RMSE (error medio cuadrático).
Las transformaciones de semejanza permiten rotar, desplazar y escalar geometrías, mientras que las afines permiten rotar, desplazar, escalar (incluso aplicando diferentes factores de escala en cada eje) y torcer geometrías. Por lo menos dos (2) puntos de control son requeridos para la de semejanza y tres (3) para la afín. Ver (Knippers, 2009) para una introducción a las fórmulas de cada transformación.
Ejemplo
Como ejemplo, podemos tomar estos dos conjuntos de datos: Departamento de La Guajira (conjunto de datos ‘A’, fuente: SIGOT) y límite de Colombia (conjunto de datos ‘B’). Los datos se incluyen solo para propósitos de demostración. Los datos comparten el sistema de referencia espacial especificado por el código EPSG 3116.
Importémoslos en R a través del paquete rgdal (primero hay que extraer los datos y ajustar la ruta a ellos cuando se requiera):
1 2 3 |
require("rgdal") l = readOGR(dsn="/path/to/the/data/", layer="La_Guajira") r = readOGR(dsn="/path/to/the/data/", layer="Colombia") |
Así es como lucen (La Guajira: línea roja, Límite de Colombia: Polígono verde con borde negro):
4 5 |
plot(l,col="red") plot(r,add=TRUE,col="lightGreen") |
Podemos también mirarlos de cerca:
6 7 |
plot(r,col="lightGreen",xlim=c(1206309,1264449),ylim=c(1832839,1864827)) plot(l,col="red",add=TRUE) |
Puede notarse que los datos simplemente no coinciden, a pesar de compartir el mismo sistema de referencia. En este ejemplo llevaremos el límite de La Guajira a coincidir con el límite de Colombia.
Puntos de control
Para reposicionar el conjunto de datos ‘A’ tomando el ‘B’ como base, seleccionamos puntos presentes en ambos conjuntos de datos y disponemos las coordenadas de esos puntos en cuatro columnas, así:
X SOURCE Y SOURCE X TARGET Y TARGET
Las coordenadas fuente (source) corresponden al conjunto de datos ‘A’, es decir, el conjunto de datos a ser transformado, mientras que las coordenadas destino (target) corresponden al conjunto de datos base (‘B’).
Los puntos de control deben estar bien distribuidos sobre la región a ser reposicionada (Iliffe, pp.135-137). Para este ejemplo, he seleccionado 16 puntos que están incluidos en el paquete ‘vec2dtransf’ como datos de ejemplo, así que podemos acceder a ellos ejecutando:
8 9 10 11 |
require("vec2dtransf") data(control.points) options(digits=10) # To see the decimals control.points |
Estos puntos de control están bien distribuidos sobre la región a ser transformada:
12 13 |
plot(control.points[2:3],pch=19) plot(l,add=TRUE,col="red") |
Transformando coordenadas
Primero definimos el objeto transformación afín a partir de los puntos de control (solo pasamos las columnas de las coordenadas):
14 |
aft = AffineTransformation(control.points[2:5]) |
Ahora calculamos los parámetros de la transformación:
15 |
calculateParameters(aft) |
Podemos ver los parámetros calculados de esta forma:
16 |
getParameters(aft) |
Y pueso que usamos más puntos de control que los requeridos para una transformación afín, podemos acceder a los residuales y al RMSE (en metros), los cuales resultan de los mínimos cuadrados.
17 18 |
getResiduals(aft) getRMSE(aft) |
Los residuales son las diferencias entre las coordenadas fuente transformadas y las coordenadas destino de los puntos de control, mientras que el RMSE mide la desviación general de las coordenadas fuente transformadas con respecto a las coordenadas destino de los puntos de control. Más información sobre residuales y RMSE puede encontrarse en (Iliffe, 2008).
Para ver el efecto global de la transformación, podríamos plotear la transformación de una grilla regular, con base en (UC Davis Soil Resource Laboratory, 2007):
19 |
plotGridTransformation(aft,bbox(l),100) |
Y agregar el conjunto de datos ‘A’:
20 |
plot(l,add=TRUE,col="red") |
Finalmente, realizamos la transformación afín sobre el conjunto de datos ‘A’:
21 |
newLines = applyTransformation(aft,l) |
Así es como lucen los datos ahora (el conjunto de datos transformado está en color azul):
22 23 |
plot(newLines,col="blue") plot(r,add=TRUE,col="lightGreen") |
Y en la vista detallada:
24 25 |
plot(r,col="lightGreen",xlim=c(1206309,1264449),ylim=c(1832839,1864827)) plot(newLines,col="blue",add=TRUE) |
Mucho mejor, ¿cierto? Comparémoslo con el conjunto de datos original (el mal posicionado):
26 |
plot(l,col="red",add=TRUE) |
Usando el paquete rgdal podemos exportar el resultado como, por ejemplo, un Shapefile, en este caso a un directorio temporal:
27 |
writeOGR(newLines, tempdir(), "transformed_Guajira", driver="ESRI Shapefile") |
Referencias
- Iliffe, J. and Lott, R. Datums and map projections: For remote sensing, GIS and surveying. Section 4.5. pp.109-117,135-137, 2008.
- UC Davis Soil Resource Laboratory. Case Study: Fixing Bad TIGER Line data with R and PostGIS. 2007. <URL: http://casoilresource.lawr.ucdavis.edu/drupal/node/433>
- Knippers, R. 2D Cartesian coordinate transformations. Section 5.4. 2009. <URL: http://kartoweb.itc.nl/geometrics/Coordinate%20transformations/coordtrans.html>
Comentarios recientes