Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save kikocorreoso/4708735 to your computer and use it in GitHub Desktop.
Save kikocorreoso/4708735 to your computer and use it in GitHub Desktop.
Usando google earth con ayuda de python y pykml (II).ipynb ipython notebook usado para la entrada http://wp.me/p2hEpj-p3 en pybonacci.wordpress.com
Display the source blob
Display the rendered blob
Raw
{
"metadata": {
"name": "Usando google earth con ayuda de python y pykml (II)"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Esta es la segunda parte de la entrada 'Usando google earth con ayuda de python y pykml (I)'. En el ejemplo de hoy vamos a ver algo m\u00e1s relacionado con la ciencia que lo visto anteriormente, que no era m\u00e1s que un aburrido ejemplo para introduciros en google earth, kml/kmz y pykml.\n",
"\n",
"**[Para esta entrada se ha usado numpy 1.6.1, matplotlib 1.1.1rc, netCDF4 1.0.2, pykml 0.1.0, lxml 2.3.2 (es un prerrequisito de pykml) en el notebook de IPython 0.13 corriendo todo sobre python 2.7.3]**\n",
"\n",
"Primero importamos una serie de cosas que usaremos."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import os\n",
"from zipfile import ZipFile\n",
"import datetime as dt\n",
"\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"import netCDF4 as nc\n",
"from lxml import etree"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En este caso vamos a representar el hurac\u00e1n Iv\u00e1n, que se gener\u00f3 en 2004 en el Atl\u00e1ntico occidental y lleg\u00f3 a convertirse en un hurac\u00e1n muy potente. Vamos a representar la evoluci\u00f3n de la presi\u00f3n al nivel del mar y las zonas de viento superiores a 50 km/h a medida que vamos viendo la trayectoria del hurac\u00e1n. El valor de viento se obtiene seg\u00fan los datos que vamos a usar, estos datos no son muy rigurosos en este sentido por lo que esto solo sirve como ejemplo.\n",
"\n",
"[Descargaremos datos](http://data-portal.ecmwf.int/data/d/interim_full_daily/) de [rean\u00e1lisis del centro europeo de predicci\u00f3n a medio plazo (ECMWF, por sus siglas en ingl\u00e9s, European Center for Medium-range Weather Forecasts)](http://es.wikipedia.org/wiki/Centro_Europeo_de_Previsiones_Meteorol%C3%B3gicas_a_Plazo_Medio) para las fechas de ocurrencia del hurac\u00e1n. Los datos usados en el ejemplo los pod\u00e9is [descargar de aqu\u00ed (al fichero le he puesto extensi\u00f3n .gif para poder subirlo a wordpress pero es un fichero netcdf, para que funcione el ejemplo deber\u00e9is quitar la extensi\u00f3n gif o cambiar el nombre del fichero en el c\u00f3digo)](http://wp.me/a2hEpj-p8).\n",
"\n",
"Ahora importamos todo lo necesario de pykml para este ejemplo:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from pykml.factory import nsmap\n",
"from pykml.factory import KML_ElementMaker as KML\n",
"from pykml.factory import GX_ElementMaker as GX"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ahora vamos a detallar la trayectoria del hurac\u00e1n con su posici\u00f3n, fecha y a en qu\u00e9 categor\u00eda se encontraba en cada fecha. La categor\u00eda del hurac\u00e1n puede ser de 1 a 5 siendo el 5 el m\u00e1s extremo, seg\u00fan la [escala Saffir-Simpson](http://www.nhc.noaa.gov/aboutsshws.php). En la categor\u00eda del hurac\u00e1n tambi\u00e9n vamos a detallar tambi\u00e9n los momentos en que no es hurac\u00e1n como tal sino que es tormenta tropical o depresi\u00f3n tropical."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Desde la versi\u00f3n 5.0, google ha a\u00f1adido una serie de extensiones\n",
"# a kml que lo hacen m\u00e1s potente e interactivo\n",
"# https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions\n",
"# Definimos una variable para el espacio de nombres de las\n",
"# extensiones de google que vamos a usar.\n",
"gxns = '{' + nsmap['gx'] + '}'\n",
"\n",
"# Ponemos la informaci\u00f3n referente al hurac\u00e1n:\n",
"\n",
"# La posici\u00f3n en longitud cada 6h desde que es tormenta tropical\n",
"lon_hur = [-30.3, -32.1, -33.6, -35, -36.5, -38.2, -39.9, \n",
" -41.4, -43.4, -45.1, -46.8, -48.5, -50.5, -52.5, \n",
" -54.4, -56.1, -57.8, -59.4, -61.1, -62.6, -64.1, \n",
" -65.5, -67, -68.3, -69.5, -70.8, -71.9, -72.8, \n",
" -73.8, -74.7, -75.8, -76.5, -77.6, -78.4, -79, \n",
" -79.6, -80.4, -81.2, -82.1, -82.8, -83.5, -84.1, \n",
" -84.7, -85.1, -85.6, -86, -86.5, -87, -87.4, \n",
" -87.9, -88.2, -88.2, -87.9, -87.7, -87.4, -86.5, \n",
" -85.7, -84, -82.3, -80.5, -78.5, -76.7, -75.5, \n",
" -74, -74, -74.5, -75.8, -77.5, -78.5, -78.7, -79.1, \n",
" -79.7, -80.6, -81.7, -82.8, -84.1, -86.1, -87.3, \n",
" -88.6, -89.5, -91, -92.2, -92.7, -93.2, -94.2]\n",
"# La posici\u00f3n en latitud cada 6h desde que es tormenta tropical\n",
"lat_hur = [9.7, 9.5, 9.3, 9.1, 8.9, 8.9, 9, 9.3, 9.5, 9.8, \n",
" 10.2, 10.6, 10.8, 11, 11.3, 11.2, 11.3, 11.6, \n",
" 11.8, 12, 12.3, 12.6, 13, 13.3, 13.7, 14.2, 14.7, \n",
" 15.2, 15.7, 16.2, 16.8, 17.3, 17.4, 17.7, 18, \n",
" 18.2, 18.4, 18.8, 19.1, 19.5, 19.9, 20.4, 20.9, \n",
" 21.6, 22.4, 23, 23.7, 24.7, 25.6, 26.7, 27.9, \n",
" 28.9, 30, 31.4, 32.5, 33.8, 34.7, 35.4, 36.2, \n",
" 37, 37.7, 38.4, 38, 37.5, 36, 34.5, 32.8, 31, \n",
" 29, 27.5, 26.4, 26.1, 25.9, 25.8, 25.2, 24.8, \n",
" 25.1, 26, 26.5, 27.1, 27.9, 28.9, 29.2, 29.6, 30.1]\n",
"# El estado del huracan cada 6h\n",
"tipo = ['TS','TS','TS','TS','TS','TS','TS',\n",
" 'TS','H1','H2','H3','H4','H3','H3',\n",
" 'H2','H2','H2','H3','H3','H4','H4',\n",
" 'H4','H4','H4','H5','H5','H4','H4',\n",
" 'H4','H4','H4','H4','H4','H4','H5',\n",
" 'H5','H4','H4','H4','H5','H5','H5',\n",
" 'H5','H5','H5','H4','H4','H4','H4',\n",
" 'H4','H4','H3','H3','H1','TS','TD',\n",
" 'TD','TD','TD','TD','TD','TD','TD',\n",
" 'TS','TS','TS','TS','TS','TS','TD',\n",
" 'TD','TD','TD','TD','TD','TD','TD',\n",
" 'TD','TD','TS','TS','TS','TS','TD','TD']\n",
"# La fecha cada 6h desde el 21/09/1998 a las 18Z\n",
"fecha_inicial = dt.datetime(2004, 9, 3, 6, 0, 0)\n",
"deltat = dt.timedelta(hours = 6)\n",
"fechas = [fecha_inicial + deltat * i for i in range(len(tipo))]\n",
"# El icono que vamos a asignar a cada 6h, en funci\u00f3n del estado de la tormenta.\n",
"iconos = {'TD':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/td.gif',\n",
" 'TS':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/ts.gif',\n",
" 'H1':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h1.gif',\n",
" 'H2':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h2.gif',\n",
" 'H3':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h3.gif',\n",
" 'H4':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h4.gif',\n",
" 'H5':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h5.gif'}"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ahora vamos a definir una funci\u00f3n que es la que usaremos para dibujar los datos de presi\u00f3n y viento que luego ir\u00e1n en el fichero kml/kmz"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Funci\u00f3n que crea las im\u00e1genes\n",
"def pinta_mapas(nombre, lon, lat, mslp, u, v, loni, lati):\n",
" fig = plt.figure()\n",
" fig.set_size_inches(8, 8. * lat.shape[0] / lon.shape[0])\n",
" ax = plt.Axes(fig, [0., 0., 1., 1.])\n",
" ax.set_axis_off()\n",
" fig.add_axes(ax)\n",
" ax.contour(lon, lat, mslp / 100., np.arange(900,1100.,2.),\n",
" colors='y',linewidths=2, aspect='normal')\n",
" wspd = np.sqrt(u * u + v * v) * 3.6\n",
" ax.contourf(lon, lat, wspd, np.arange(50, 250, 25), \n",
" colors = plt.hot(), aspect='normal')\n",
" plt.savefig(nombre, dpi = 80, transparent=True)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Y ahora vamos a leer los datos que os hay\u00e1is descargado:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Leemos los datos de presi\u00f3n, velocidad, latitud y longitud.\n",
"data = nc.Dataset('Usando google earth con ayuda de python y pykml (II)/data.nc')\n",
"lon = data.variables['longitude'][:]\n",
"lat = data.variables['latitude'][:]\n",
"# presi\u00f3n = mslp por mean sea level pressure,\n",
"# presi\u00f3n media al nivel del mar\n",
"mslp = data.variables['msl'][1:-2,:,:]\n",
"lon[lon > 180] = lon - 360\n",
"u = data.variables['u10'][1:-2,:,:]\n",
"v = data.variables['v10'][1:-2,:,:] "
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Y ahora es cuando vamos a crear el fichero kml haciendo un bucle donde se van a ir a\u00f1adiendo cosas tanto a la carpeta (folder) como al Tour. Voy explic\u00e1ndo lo que hace cada cosa mediante comentarios en el c\u00f3digo:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Como vamos a crear un fichero kmz que incluir\u00e1 todas las im\u00e1genes\n",
"# abrimos el fichero kmz donde iremos guard\u00e1ndolo todo\n",
"fich_kmz = ZipFile('HuracanIvan.kmz', 'w')\n",
"\n",
"# Desde la versi\u00f3n 5.0, google ha a\u00f1adido una serie de extensiones\n",
"# a kml que lo hacen m\u00e1s potente e interactivo\n",
"# https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions\n",
"# Definimos una variable para el espacio de nombres de las\n",
"# extensiones de google que vamos a usar.\n",
"gxns = '{' + nsmap['gx'] + '}'\n",
"# start with a base KML tour and playlist\n",
"fich_kml = KML.kml(\n",
" KML.Document(\n",
" GX.Tour(\n",
" KML.name(u\"\u00a1Reprod\u00faceme!\"),\n",
" GX.Playlist(),\n",
" ),\n",
" KML.Folder(\n",
" KML.name('Ruta del huracan Ivan'),\n",
" id='lugares',\n",
" ),\n",
" )\n",
")\n",
"# Hacemos un bucle para recorrer todos los 'momentos'\n",
"# del hurac\u00e1n definidos en el bloque de c\u00f3digo anterior\n",
"for i in range(len(tipo)):\n",
" # Primero volamos hasta cada posici\u00f3n del hurac\u00e1n (cada 6h)\n",
" # y nos quedamos observ\u00e1ndolo desde el espacio \n",
" # A unos 2000 km de altura. \n",
" # Para ello usamos FlyTo donde especificamos\n",
" # La duraci\u00f3n en segundos, si queremos que vaya m\u00e1s r\u00e1pido o m\u00e1s\n",
" # lento, el modo de vuelo, que puede ser m\u00e1s suave (smooth) o \n",
" # m\u00e1s r\u00e1pido (bounce) y la posici\u00f3n.\n",
" # La posici\u00f3n la determinan longitude, latitude y altitude,\n",
" # tilt es el \u00e1ngulo desde la vertical (en este caso lo vemos\n",
" # verticalmente (0\u00ba)) y range es la distancia desde la posici\u00f3n\n",
" # determinada por longitude, latitude y altitude, teniendo en\n",
" # cuenta el \u00e1ngulo tilt (en este caso hemos puesto 2.000 km)\n",
" # Las unidades para latitude y longitude son grados, para altitude\n",
" # y range son metros y para tilt son grados desde la vertical.\n",
" # En el gr\u00e1fico de este enlace se ver\u00e1 mejor:\n",
" # https://developers.google.com/kml/documentation/images/lookAt.gif\n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(\n",
" GX.FlyTo(\n",
" GX.duration(1),\n",
" GX.flyToMode(\"bounce\"),\n",
" KML.LookAt(\n",
" KML.longitude(lon_hur[i]),\n",
" KML.latitude(lat_hur[i]),\n",
" KML.altitude(0),\n",
" KML.heading(0),\n",
" KML.tilt(0),\n",
" KML.range(10000000.0),\n",
" KML.altitudeMode(\"relativeToGround\"),\n",
" )\n",
" ),\n",
" )\n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(GX.Wait(GX.duration(0.1)))\n",
" # A\u00f1adimos informaci\u00f3n de cada lugar en la carpeta\n",
" # como 'placemarks' o marcas de posici\u00f3n. En el nombre y en la descripci\u00f3n\n",
" # del placemark se puede meter HTML y CSS, por si quer\u00e9is ser un poco m\u00e1s\n",
" # creativos y dejar la informaci\u00f3n de una forma visible y bonita.\n",
" # extrude == 1 indica que el nombre del lugar se har\u00e1 visible. Si queremos\n",
" # hacerlo visible solo cuando lo seleccionemos pues le ponemos el valor 0\n",
" desc = '<![CDATA[{0}<br>lon: {1}<br>lat: {2}<br>Estado del huracan: {3}]]>'\n",
" desc = desc.format(fechas[i].isoformat()[0:13], lon_hur[i], lat_hur[i], tipo[i])\n",
" fich_kml.Document.Folder.append(\n",
" KML.Placemark(\n",
" KML.name(fechas[i].isoformat()[0:13]),\n",
" KML.description(\n",
" \"<h1>{0}</h1><br/>{1}\".format('IVAN', desc,)\n",
" ),\n",
" KML.Style(\n",
" KML.IconStyle(\n",
" KML.Icon(KML.href(iconos[tipo[i]]),)\n",
" ),\n",
" ),\n",
" KML.Point(\n",
" KML.extrude(1),\n",
" KML.altitudeMode(\"relativeToGround\"),\n",
" KML.coordinates(\"{0},{1},0\".format(lon_hur[i],\n",
" lat_hur[i],)\n",
" ),\n",
" ),\n",
" id=fechas[i].isoformat()[0:13]\n",
" )\n",
" )\n",
" # Aqu\u00ed le indicamos a la reproducci\u00f3n que nos \n",
" # muestre el 'globo' o 'vi\u00f1eta' con la informaci\u00f3n del sitio,\n",
" # nombre y descripci\u00f3n. Visibility a 1 para \n",
" # hacerlo visible.\n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(\n",
" GX.AnimatedUpdate(\n",
" GX.duration(1.0),\n",
" KML.Update(\n",
" KML.targetHref(),\n",
" KML.Change(\n",
" KML.Placemark(\n",
" KML.visibility(1),\n",
" GX.balloonVisibility(1),\n",
" targetId=fechas[i].isoformat()[0:13]\n",
" )\n",
" )\n",
" )\n",
" )\n",
" )\n",
" # creamos el mapa\n",
" nombre = '{0}.png'.format(fechas[i].isoformat()[0:13])\n",
" pinta_mapas(nombre, lon, lat, mslp[i,:,:], u[i,:,:], v[i,:,:], lon_hur[i], lat_hur[i])\n",
" # a\u00f1adimos el mapa al fichero kmz\n",
" fich_kmz.write(nombre)\n",
" # y eliminamos el fichero png\n",
" os.remove(nombre)\n",
" # a\u00f1adimos la ruta de las im\u00e1genes en el \n",
" # fichero kml (que estar\u00e1 contenido dentro\n",
" # del fichero kmz.\n",
" fich_kml.Document.Folder.append(\n",
" KML.GroundOverlay(\n",
" KML.name(nombre),\n",
" KML.visibility(0),\n",
" KML.Icon(KML.href(nombre)),\n",
" KML.LatLonBox(KML.north(np.max(lat)),\n",
" KML.south(np.min(lat)),\n",
" KML.east(np.max(lon)),\n",
" KML.west(np.min(lon)),\n",
" ),\n",
" id = nombre\n",
" )\n",
" ) \n",
" \n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(GX.Wait(GX.duration(0.1)))\n",
" # Aqu\u00ed le indicamos a la reproducci\u00f3n que nos \n",
" # muestre el 'mapa'. Visibility a 1 para \n",
" # hacerlo visible.\n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(\n",
" GX.AnimatedUpdate(\n",
" GX.duration(1.0),\n",
" KML.Update(\n",
" KML.targetHref(),\n",
" KML.Change(\n",
" KML.GroundOverlay(\n",
" KML.visibility(1),\n",
" #GX.balloonVisibility(1),\n",
" targetId=nombre\n",
" )\n",
" )\n",
" )\n",
" )\n",
" )\n",
" # Esta parte de c\u00f3digo har\u00e1 que desaparezca el 'globo'\n",
" # con la informaci\u00f3n as\u00ed no tendremos el 'globo' \n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(\n",
" GX.AnimatedUpdate(\n",
" GX.duration(0.1),\n",
" KML.Update(\n",
" KML.targetHref(),\n",
" KML.Change(\n",
" KML.Placemark(\n",
" GX.balloonVisibility(0),\n",
" targetId=fechas[i].isoformat()[0:13]\n",
" )\n",
" )\n",
" )\n",
" )\n",
" )\n",
" # Esta parte de c\u00f3digo har\u00e1 que desaparezca el 'mapa'\n",
" fich_kml.Document[gxns+\"Tour\"].Playlist.append(\n",
" GX.AnimatedUpdate(\n",
" GX.duration(1.0),\n",
" KML.Update(\n",
" KML.targetHref(),\n",
" KML.Change(\n",
" KML.GroundOverlay(\n",
" KML.visibility(0),\n",
" targetId=nombre\n",
" )\n",
" )\n",
" )\n",
" )\n",
" )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"En el siguiente c\u00f3digo lo que hacemos es guardar el fichero kml final. Si os acord\u00e1is, en alg\u00fan momento hemos puesto algunas etiquetas HTML. Algunos s\u00edmbolos especiales se han transformado por lo que los volvemos a dejar como estaban ya que si no google earth no entender\u00e1 ese kml correctamente."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# guardamos toda la informaci\u00f3n\n",
"# del kml en un fichero kml\n",
"outfile = file('HuracanIvan.kml','w')\n",
"salida = etree.tostring(fich_kml, pretty_print=True)\n",
"salida = salida.replace('&lt;', '<')\n",
"salida = salida.replace('&gt;', '>')\n",
"outfile.write(salida)\n",
"outfile.close()\n",
"# guardamos el fichero kml dentro del\n",
"# fichero kmz\n",
"fich_kmz.write('HuracanIvan.kml')\n",
"# Eliminamos el fichero kml\n",
"os.remove('HuracanIvan.kml')\n",
"# Cerramos el fichero kmz\n",
"fich_kmz.close()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Por \u00faltimo, lanzamos google earth con el fichero que acabamos de crear. Las siguientes tres l\u00edneas funcionan en linux. Ahora mismo no tengo un windows cerca pero creo que deber\u00e9is cambiar `googleearth` por `googleearth.exe` y tener en cuenta los '\\' en lugar de '/' en la ruta al fichero."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"ruta = os.getcwd()\n",
"rutafich = ruta + '/HuracanIvan.kmz'\n",
"os.system('googleearth {0}'.format(rutafich))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A la izquierda de la pantalla de google earth, en mis 'lugares', dentro de 'lugares temporales' podr\u00e9is ver 'HuracanIvan.kmz'. Si lo despleg\u00e1is ver\u00e9is 'Reprod\u00faceme'. Si hac\u00e9is doble click sobre 'Reprod\u00faceme' podr\u00e9is ver algo parecido a lo siguiente:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from IPython.display import YouTubeVideo\n",
"YouTubeVideo('TPKZo4ez1V0')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
" <iframe\n",
" width=\"400\"\n",
" height=\"300\"\n",
" src=\"http://www.youtube.com/embed/TPKZo4ez1V0\"\n",
" frameborder=\"0\"\n",
" allowfullscreen\n",
" ></iframe>\n",
" "
],
"output_type": "pyout",
"prompt_number": 1,
"text": [
"<IPython.lib.display.YouTubeVideo at 0x2900710>"
]
}
],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"El resultado tampoco es maravilloso pero la idea no era el resultado, la idea era mostrar el uso de pykml e introducir un poco de kml/kmz a los que les pueda resultar \u00fatil.\n",
"\n",
"Por \u00faltimo, como coment\u00e9 en el anterior cap\u00edtulo, os dejo un enlace con un uso creativo de google earth [1].\n",
"\n",
"[1] - [http://www.fosslc.org/drupal/content/pykml-python-kml-library](http://www.fosslc.org/drupal/content/pykml-python-kml-library)"
]
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment