Traspaso de imágenes de GMT a ArcGis

Para poder traspasar información desde distintos programas, se necesita conocer dos cosas. Primero conocer bien los distintos tipos de formatos con que trabajan los programas, y segundo saber exportar los datos con las herramientas que los mismos programas proporcionan. En una próxima publicación voy a ondar sobre los distintos de formatos(grd, xyz, tiff, etc).

En esta entrada abordaremos solo el traspaso de imágenes hechas con GMT a Arcgis. La idea es crear algo, que al ser ingresado a arcgis, quede automáticamente georeferenciado. Para esto la imagen creada por GMT debe incorporar los datos dentro.

Primero hay que crear una imagen cualquiera con GMT, yo voy a hacer un sencillo mapa topografico de la zona "-70.69/-70.46/-28.58/-28.42"

Primero declaro las variables:    

#! /bin/bash#Variables

zona="-70.69/-70.46/-28.58/-28.42"grdfile="campo.grd"
proyeccion="M13c"
paleta="DEM_screen.cpt"
BE="a0.025f0.025g0.025SENW"azimut="270"
psfile="topo.ps"
#Aquí creamos el archivo para la iluminación y ajustamos la paleta a nuestra topografía
#creacion de archivos temporales

grdgradient $grdfile -Giluminacion.int -A$azimut -Nt -M grd2cpt $grdfile -C$paleta -Z > paleta2.cpt
paleta2="paleta2.cpt"

#comandos
grdimage $grdfile -Bg2m --BASEMAP_TYPE=inside -C$paleta2 -J$proyeccion -Iiluminacion.int -R$zona -V -K > $psfile
grdcontour $grdfile -J$proyeccion -C200 -A200+k0/0/0 -R$zona -W0.1p/0/0/0 -Gd10c -O >> $psfile

#Borrar archivos temporales

rm $paleta2 iluminacion.int

Es importante notar --BASEMAP_TYPE=inside, pues si agregamos un marco al mapa, al ingresarlo a ArcGis este entenderá que el marco esta dentro de la imagen, distorsionando el mapa. Luego ubicados en la carpeta donde se encuentre nuestro archivo .ps, ejecutamos el siguiente comando.

ps2raster cotas.ps -Tt -Wg -V

Este comando transformara nuestra imagen a un raster,

• -Tt: indica el formato de salida, en este caso la salida es .tiff
• -Wg: es un poco más complicado, pero en palabras simples indica que el archivo .tiff, pueda ser leido por otros programas como un "geotiff"
• -V: verbose

Para más datos sobre ps2raster, se puede teclear man ps2raster en la terminal. Es importante notar que el GMT crea un archivo .tfw, es importante que conserven este archivo y lo guarden en la misma carpeta que el archivo .tiff creado. 

Una vez creado tenemos que ingresarlo a ArcGis, yo como uso Debian GNU/Linux, corro ArcGis en una maquina virtual. 

El problema al ingresarlo es que aunque puede ser leido por ArcGis, este no reconoce el sistema cartografico inmediatamente. Lo que hice yo fue ir probando las diferentes opciones que te da arcgis hasta que me resulto. Para que no se me olvidara guarde el sistema de georeferencia como un sistema de coordenadas personalizado. 

Les dejo aquí un link para descargarlo.

Proyección GMT

Ahora solo cuando seleccionen el tipo de proyección, buscan este archivo.

Eso es todo por hoy. 

Páginas para descargar topografia y datos de sismos

Para hacer mapas topograficos en GMT yo obtengo la información de dos páginas, GEBCO08:

https://www.bodc.ac.uk/data/online_delivery/gebco/

y de ASTER GDEM


http://www.gdem.aster.ersdac.or.jp/search.jsp


La primera tiene una resolución de 30 segundos y la segunda una resolucion de 30 metros, por lo que uso la primera para mapas de zonas extensas y la segunda para zonas locales. 


Para plotear los datos de sismos uso el comando psmeca, por lo que es comodo usar globalCMT pues te entrega los datos directamente en este formato.


http://www.globalcmt.org/CMTsearch.html

Tipos de proyecciones

En GMT la opción -J establece que tipo de proyección usaremos para graficar nuestro mapa, si utilizamos el comando man psbasemap y vamos donde se especifica la opción -J se puede apreciar todas las proyecciones que permite usar GMT. Si lo han hecho o lo acaban de hacer se puede dar cuenta la cantidad de opciones que tiene.
Para intentar aclarar esto me dedicare en este post a explicar los principales tipos de proyecciones que existen y las caracteristicas de cada una.

Primero que nada, las proyecciones pueden ser equivalentes o conformes.:
Las proyecciones equivalentes son aquellas que mantienen constante la superficie a lo largo y ancho del mapa. Sin embargo las formas son distorcionadas. Ejemplo de estas proyecciones son la proyección de Albers y Mollweide.
Por otro lado las proyecciones conformales se caracterizan por mantener las formas pero pierden en el tamaño de la superficie. Ejemplos de estas proyecciones son la proyección de Mercator y Miller.

Además según como se proyectan se tienen 4 tipos principales:

Cilindricas: son proyecciones conformales, que se obtienen al proyectar la forma de la tierra en un cilindro para luego estirarlo. Se utliza sobre todo en mapamundis La distorción aumenta hacia los polos hasta el punto en que, por ejemplo, Africa y Groenlandia se ven del mismo tamaño cuando en realidad africa es más o menos 10 veces más grande que Groenlandia.Son proyecciones cilíndricas: Mercator, Lambert

Elipticas: son proyecciones con forma de elipse. Presentan los paralelos rectos y los meridianos curvos. El punto central es el punto de no-distorción y esta aumenta hacia los polos.

Cónicas: Se proyecta la superficie sobre uno o más conos que estan ubicado tangente a la superficie de la tierra. Generalmente se usa esta proyección para ver los paises y continentes.

Acimutales: Se proyecta la tierra sobre un plano que es tangente a un punto del globo. La distorción aumenta al alejarse del punto tangencial. En este tipo de proyección solo es posible mostrar un plano.

La proxima vez describire las opciones de -J que trae GMT.

Referencias:
Jorge Fallas. 2003. Proyecciones cartográficas y datum ¿Que son y para que sirven. Laboratorio de teledetección y sistemas de información geografica.
http://es.wikipedia.org/wiki/Proyecci%C3%B3n_cartogr%C3%A1fica

Introducción y animación con GMT

La verdad me es difícil señalar para que tipo de lector esta definido este blog. Tanto en conocimientos previos como la utilidad que le pueda sonsacar en su lectura. Sin embargo el blog avanzara igual, para nadie o para todos, dirigido a mi lector imaginario o quizá para un lector fugaz que encuentre aquí algo que le pueda ser de utilidad. 

No explicare GMT desde 0 pues me sería demasiado tedioso hacer esto. Por otro lado yo mismo llevo trabajando con él solo durante un par de meses así que cada vez que descubra una nueva funcionalidad de este lo escribiré en este blog. Hoy vamos a ver como a partir de una serie de imagenes 3D, en este caso topografía, vamos a conseguir hacer un sencillo vídeo(aunque de alta resolución).    

Primero vamos a obtener nuestro archivo de topografía en el que vamos a trabajar, para esto el mejor sitio que conozco para descargar (de forma gratuita) es http://www.gdem.aster.ersdac.or.jp/search.jsp. No explicare como descargarlo pues la interfaz de la pagina me parece bastante intuitiva(de todas formas si alguien tiene alguna duda que me envíe un mail o algo).  La resolución de la topografía es de 30 m. 

De aquí se pueden obtener archivos .tiff los cuales no nos sirven para tratarlos con GMT. Para solucionar esto usamos el siguiente comando en la terminal(utilizamos el archivo que dice dem):

gdal_translate -of gmt archivo.tf archivo.grd ; para más información de gdal, poner man gdal_translate en la terminal.

Ahora que tenemos nuestro archivo .grd vamos a realizar un script para diseñar nuestra imagen en 3D. No me detendré mucho en esta parte, solo mostrare los comandos necesarios. 

Primero preparamos el archivo, definimos la iluminación de la imagen 3D y luego re definimos los limites de la paleta de colores según nuestra topografía :

grdgradient archivo.grd -Ggradiente.int -A250 -Nt -M   

grd2cpt archivo.grd -Cpaleta -Z  > paleta2.cpt ; se pueden descargar paletas desde 

http://soliton.vm.bytemark.co.uk/pub/cpt-city/

Ya que este todo listo, en el script antes de poner grdview para crear nuestra imagen utilizamos el comando for para crear un ciclo:

for x in {0..360} 

do

Con esto básicamente le decimos  a la maquina que realice 361 veces los siguientes comandos para la zona longitud -70.73/-70.53, latitud -28.5/-28.33 : 

grdview archivo.grd -JM13c -JZ1.16c -Cpaleta2.cpt -R-70.73/-70.53/-28.5/-28.33 -E$x/35+v8/13 -Igradiente.int -B0.05a0.05 -Qs -P -V -K  > archivo.ps

psscale -D15c/4c/5c/0.5c -Cpaleta2.cpt -Ba400f200/a400f200:"(m)": -O >> archivo.ps

gs -sDEVICE=jpeg -dJPEGQ=100 -dNOPAUSE -dBATCH -dSAFER -r300 -sOutputFile="$x".jpg archivo.ps


rm archivo.ps


done

De esta forma mientras x varia de 0 a 360 en la parte -E$x/35+v8/13, el angulo en el cual vemos la imagen 3D ira variando, por lo que con eso se generaran 361 imágenes en la cuál la imagen ira rotando con el centro ubicado por +v8/13.

El penúltimo comando sirve para transformar los archivos .ps que crea GMT a .jpg. El ultimo sirve para eliminar los archivos .ps que ya no usaremos.

Por último luego de tener todas las imágenes las vamos a juntar con 

mencoder "mf://*.jpg" -mf fps=25 -o archivovideo.avi -ovc lavc -lavcopts vcodec=mpeg4

Para lograr un vídeo fluido recomiendo tener sobre 25 imágenes por segundo. Aquí puse 25 en donde dice fps=25. Para más detalles ver en 

http://www.mplayerhq.hu/DOCS/HTML/es/menc-feat-enc-images.html

De esta forma desarrollamos un vídeo como el siguiente :  

Utilizando este método con ciclos para crear multi-imágenes y luego unirlas en un vídeo podemos crear un sin fin de vídeos con GMT.  Podemos cambiar la elevación, la localidad, colores, etc para lograr un bonito vídeo. Espero que le sirva a alguien. Cualquier duda pueden contactarse conmigo.