import os from zipfile import ZipFile import datetime as dt import numpy as np from matplotlib import pyplot as plt import netCDF4 as nc from lxml import etree from pykml.factory import nsmap from pykml.factory import KML_ElementMaker as KML from pykml.factory import GX_ElementMaker as GX # Desde la versión 5.0, google ha añadido una serie de extensiones # a kml que lo hacen más potente e interactivo # https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions # Definimos una variable para el espacio de nombres de las # extensiones de google que vamos a usar. gxns = '{' + nsmap['gx'] + '}' # Ponemos la información referente al huracán: # La posición en longitud cada 6h desde que es tormenta tropical lon_hur = [-30.3, -32.1, -33.6, -35, -36.5, -38.2, -39.9, -41.4, -43.4, -45.1, -46.8, -48.5, -50.5, -52.5, -54.4, -56.1, -57.8, -59.4, -61.1, -62.6, -64.1, -65.5, -67, -68.3, -69.5, -70.8, -71.9, -72.8, -73.8, -74.7, -75.8, -76.5, -77.6, -78.4, -79, -79.6, -80.4, -81.2, -82.1, -82.8, -83.5, -84.1, -84.7, -85.1, -85.6, -86, -86.5, -87, -87.4, -87.9, -88.2, -88.2, -87.9, -87.7, -87.4, -86.5, -85.7, -84, -82.3, -80.5, -78.5, -76.7, -75.5, -74, -74, -74.5, -75.8, -77.5, -78.5, -78.7, -79.1, -79.7, -80.6, -81.7, -82.8, -84.1, -86.1, -87.3, -88.6, -89.5, -91, -92.2, -92.7, -93.2, -94.2] # La posición en latitud cada 6h desde que es tormenta tropical lat_hur = [9.7, 9.5, 9.3, 9.1, 8.9, 8.9, 9, 9.3, 9.5, 9.8, 10.2, 10.6, 10.8, 11, 11.3, 11.2, 11.3, 11.6, 11.8, 12, 12.3, 12.6, 13, 13.3, 13.7, 14.2, 14.7, 15.2, 15.7, 16.2, 16.8, 17.3, 17.4, 17.7, 18, 18.2, 18.4, 18.8, 19.1, 19.5, 19.9, 20.4, 20.9, 21.6, 22.4, 23, 23.7, 24.7, 25.6, 26.7, 27.9, 28.9, 30, 31.4, 32.5, 33.8, 34.7, 35.4, 36.2, 37, 37.7, 38.4, 38, 37.5, 36, 34.5, 32.8, 31, 29, 27.5, 26.4, 26.1, 25.9, 25.8, 25.2, 24.8, 25.1, 26, 26.5, 27.1, 27.9, 28.9, 29.2, 29.6, 30.1] # El estado del huracan cada 6h tipo = ['TS','TS','TS','TS','TS','TS','TS', 'TS','H1','H2','H3','H4','H3','H3', 'H2','H2','H2','H3','H3','H4','H4', 'H4','H4','H4','H5','H5','H4','H4', 'H4','H4','H4','H4','H4','H4','H5', 'H5','H4','H4','H4','H5','H5','H5', 'H5','H5','H5','H4','H4','H4','H4', 'H4','H4','H3','H3','H1','TS','TD', 'TD','TD','TD','TD','TD','TD','TD', 'TS','TS','TS','TS','TS','TS','TD', 'TD','TD','TD','TD','TD','TD','TD', 'TD','TD','TS','TS','TS','TS','TD','TD'] # La fecha cada 6h desde el 21/09/1998 a las 18Z fecha_inicial = dt.datetime(2004, 9, 3, 6, 0, 0) deltat = dt.timedelta(hours = 6) fechas = [fecha_inicial + deltat * i for i in range(len(tipo))] # El icono que vamos a asignar a cada 6h, en función del estado de la tormenta. iconos = {'TD':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/td.gif', 'TS':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/ts.gif', 'H1':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h1.gif', 'H2':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h2.gif', 'H3':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h3.gif', 'H4':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h4.gif', 'H5':'http://www.srh.noaa.gov/gis/kml/hurricanetrack/hurrimages/h5.gif'} # Función que crea las imágenes def pinta_mapas(nombre, lon, lat, mslp, u, v, loni, lati): fig = plt.figure() fig.set_size_inches(8, 8. * lat.shape[0] / lon.shape[0]) ax = plt.Axes(fig, [0., 0., 1., 1.]) ax.set_axis_off() fig.add_axes(ax) ax.contour(lon, lat, mslp / 100., np.arange(900,1100.,2.), colors='y',linewidths=2, aspect='normal') wspd = np.sqrt(u * u + v * v) * 3.6 ax.contourf(lon, lat, wspd, np.arange(50, 250, 25), colors = plt.hot(), aspect='normal') plt.savefig(nombre, dpi = 80, transparent=True) # Leemos los datos de presión, velocidad, latitud y longitud. data = nc.Dataset('Usando google earth con ayuda de python y pykml (II)/data.nc') lon = data.variables['longitude'][:] lat = data.variables['latitude'][:] # presión = mslp por mean sea level pressure, # presión media al nivel del mar mslp = data.variables['msl'][1:-2,:,:] lon[lon > 180] = lon - 360 u = data.variables['u10'][1:-2,:,:] v = data.variables['v10'][1:-2,:,:] # Como vamos a crear un fichero kmz que incluirá todas las imágenes # abrimos el fichero kmz donde iremos guardándolo todo fich_kmz = ZipFile('HuracanIvan.kmz', 'w') # Desde la versión 5.0, google ha añadido una serie de extensiones # a kml que lo hacen más potente e interactivo # https://developers.google.com/kml/documentation/kmlreference?hl=es#kmlextensions # Definimos una variable para el espacio de nombres de las # extensiones de google que vamos a usar. gxns = '{' + nsmap['gx'] + '}' # start with a base KML tour and playlist fich_kml = KML.kml( KML.Document( GX.Tour( KML.name(u"¡Reprodúceme!"), GX.Playlist(), ), KML.Folder( KML.name('Ruta del huracan Ivan'), id='lugares', ), ) ) # Hacemos un bucle para recorrer todos los 'momentos' # del huracán definidos en el bloque de código anterior for i in range(len(tipo)): # Primero volamos hasta cada posición del huracán (cada 6h) # y nos quedamos observándolo desde el espacio # A unos 2000 km de altura. # Para ello usamos FlyTo donde especificamos # La duración en segundos, si queremos que vaya más rápido o más # lento, el modo de vuelo, que puede ser más suave (smooth) o # más rápido (bounce) y la posición. # La posición la determinan longitude, latitude y altitude, # tilt es el ángulo desde la vertical (en este caso lo vemos # verticalmente (0º)) y range es la distancia desde la posición # determinada por longitude, latitude y altitude, teniendo en # cuenta el ángulo tilt (en este caso hemos puesto 2.000 km) # Las unidades para latitude y longitude son grados, para altitude # y range son metros y para tilt son grados desde la vertical. # En el gráfico de este enlace se verá mejor: # https://developers.google.com/kml/documentation/images/lookAt.gif fich_kml.Document[gxns+"Tour"].Playlist.append( GX.FlyTo( GX.duration(1), GX.flyToMode("bounce"), KML.LookAt( KML.longitude(lon_hur[i]), KML.latitude(lat_hur[i]), KML.altitude(0), KML.heading(0), KML.tilt(0), KML.range(10000000.0), KML.altitudeMode("relativeToGround"), ) ), ) fich_kml.Document[gxns+"Tour"].Playlist.append(GX.Wait(GX.duration(0.1))) # Añadimos información de cada lugar en la carpeta # como 'placemarks' o marcas de posición. En el nombre y en la descripción # del placemark se puede meter HTML y CSS, por si queréis ser un poco más # creativos y dejar la información de una forma visible y bonita. # extrude == 1 indica que el nombre del lugar se hará visible. Si queremos # hacerlo visible solo cuando lo seleccionemos pues le ponemos el valor 0 desc = 'lon: {1}
lat: {2}
Estado del huracan: {3}]]>' desc = desc.format(fechas[i].isoformat()[0:13], lon_hur[i], lat_hur[i], tipo[i]) fich_kml.Document.Folder.append( KML.Placemark( KML.name(fechas[i].isoformat()[0:13]), KML.description( "

{0}


{1}".format('IVAN', desc,) ), KML.Style( KML.IconStyle( KML.Icon(KML.href(iconos[tipo[i]]),) ), ), KML.Point( KML.extrude(1), KML.altitudeMode("relativeToGround"), KML.coordinates("{0},{1},0".format(lon_hur[i], lat_hur[i],) ), ), id=fechas[i].isoformat()[0:13] ) ) # Aquí le indicamos a la reproducción que nos # muestre el 'globo' o 'viñeta' con la información del sitio, # nombre y descripción. Visibility a 1 para # hacerlo visible. fich_kml.Document[gxns+"Tour"].Playlist.append( GX.AnimatedUpdate( GX.duration(1.0), KML.Update( KML.targetHref(), KML.Change( KML.Placemark( KML.visibility(1), GX.balloonVisibility(1), targetId=fechas[i].isoformat()[0:13] ) ) ) ) ) # creamos el mapa nombre = '{0}.png'.format(fechas[i].isoformat()[0:13]) pinta_mapas(nombre, lon, lat, mslp[i,:,:], u[i,:,:], v[i,:,:], lon_hur[i], lat_hur[i]) # añadimos el mapa al fichero kmz fich_kmz.write(nombre) # y eliminamos el fichero png os.remove(nombre) # añadimos la ruta de las imágenes en el # fichero kml (que estará contenido dentro # del fichero kmz. fich_kml.Document.Folder.append( KML.GroundOverlay( KML.name(nombre), KML.visibility(0), KML.Icon(KML.href(nombre)), KML.LatLonBox(KML.north(np.max(lat)), KML.south(np.min(lat)), KML.east(np.max(lon)), KML.west(np.min(lon)), ), id = nombre ) ) fich_kml.Document[gxns+"Tour"].Playlist.append(GX.Wait(GX.duration(0.1))) # Aquí le indicamos a la reproducción que nos # muestre el 'mapa'. Visibility a 1 para # hacerlo visible. fich_kml.Document[gxns+"Tour"].Playlist.append( GX.AnimatedUpdate( GX.duration(1.0), KML.Update( KML.targetHref(), KML.Change( KML.GroundOverlay( KML.visibility(1), #GX.balloonVisibility(1), targetId=nombre ) ) ) ) ) # Esta parte de código hará que desaparezca el 'globo' # con la información así no tendremos el 'globo' fich_kml.Document[gxns+"Tour"].Playlist.append( GX.AnimatedUpdate( GX.duration(0.1), KML.Update( KML.targetHref(), KML.Change( KML.Placemark( GX.balloonVisibility(0), targetId=fechas[i].isoformat()[0:13] ) ) ) ) ) # Esta parte de código hará que desaparezca el 'mapa' fich_kml.Document[gxns+"Tour"].Playlist.append( GX.AnimatedUpdate( GX.duration(1.0), KML.Update( KML.targetHref(), KML.Change( KML.GroundOverlay( KML.visibility(0), targetId=nombre ) ) ) ) ) # guardamos toda la información # del kml en un fichero kml outfile = file('HuracanIvan.kml','w') salida = etree.tostring(fich_kml, pretty_print=True) salida = salida.replace('<', '<') salida = salida.replace('>', '>') outfile.write(salida) outfile.close() # guardamos el fichero kml dentro del # fichero kmz fich_kmz.write('HuracanIvan.kml') # Eliminamos el fichero kml os.remove('HuracanIvan.kml') # Cerramos el fichero kmz fich_kmz.close() ruta = os.getcwd() rutafich = ruta + '/HuracanIvan.kmz' os.system('googleearth {0}'.format(rutafich)) from IPython.display import YouTubeVideo YouTubeVideo('TPKZo4ez1V0')