Continuando lo que enseñó Juanlu en la anterior entrada vamos a mostrar líneas de nivel y temperatura del aire en la superficie, en este caso la presión al nivel del mar del día 01 de enero de 2012 a las 00.00 UTC según los datos extraídos del reanálisis NCEP/NCAR, sobre un mapa con la ayuda de la librería Basemap.
Como los datos del reanálisis NCEP/NCAR vienen en formato netCDF usaremos la librería netcdf4-python. El formato netCDF es un estándar abierto y es ampliamente usado en temas de ciencias de la tierra, atmósfera, climatología, meteorología,… No es estrictamente necesario usar netcdf4-python para acceder a ficheros netCDF puesto que desde scipy tenéis esta funcionalidad. Pero bueno, yo uso esta por una serie de ventajas que veremos otro día.
En la presente entrada se ha usado python 2.7.2, numpy 1.6.1, matplotlib 1.1.0, netCDF4 0.9.7 y Basemap 1.0.2.
Primero de todo vamos a importar todo lo que necesitamos:
[sourcecode language=”python”]
## Importamos las librerías que nos hacen falta
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits import basemap as bm
[/sourcecode]
Los ficheros netCDF de presión al nivel del mar y de Temperatura del aire de la superficie los podéis descargar de aquí y aquí, respectivamente. Veréis un enlace que pone ‘FTP a copy of the file’, lo pincháis y guardáis en el mismo sitio donde tengáis el script que estamos haciendo en la presente entrada.
Una vez que tenemos los ficheros los podemos abrir usando la librería netCDF4-python:
[sourcecode language=”python”]
## Abrimos los ficheros de datos,
## el nombre de los ficheros lo tendréis que cambiar
## con el nombre de los ficheros que os habéis descargado
slp = nc.Dataset(‘X83.34.8.250.104.4.18.19.nc’) #slp por ‘sea level pressure’
tsfc = nc.Dataset(‘X83.34.8.250.104.4.15.31.nc’) #tsfc ‘por temperature at surface’
[/sourcecode]
1 |
Para saber las variables que tenemos en cada fichero netCDF podemos escribir lo siguiente:
[sourcecode language=”python”]
## Qué variables hay dentro de cada netCDF
print slp.variables
print tsfc.variables
[/sourcecode]
El output que veremos para slp será:
[sourcecode language=”python”]
OrderedDict([(u’lat’, <netCDF4.Variable object at 0x05F72FA8>), (u’lon’, <netCDF4.Variable object at 0x0603C198>), (u’time’, <netCDF4.Variable object at 0x0603C150>), (u’slp’, <netCDF4.Variable object at 0x060A1FA8>)])
[/sourcecode]
De la misma forma, para tsfc tendremos:
[sourcecode language=”python”]
OrderedDict([(u’lat’, <netCDF4.Variable object at 0x060AE108>), (u’lon’, <netCDF4.Variable object at 0x060AE030>), (u’time’, <netCDF4.Variable object at 0x060AE078>), (u’air’, <netCDF4.Variable object at 0x060AE150>)])
[/sourcecode]
Nos interesa acceder a las variables ‘slp’, ‘air’, ‘lat’ y ‘lon’. Las dos últimas variables son las mismas tanto en slp como en tsfc ya que nos hemos descargado el campo de presión y temperatura para la misma fecha y para la misma área. Por tanto, para acceder a los datos y poder dibujarlos hacemos lo siguiente:
[sourcecode language=”python”]
slpdata = slp.variables[‘slp’][:] #Obtenemos los datos en Pa
tsfcdata = tsfc.variables[‘air’][:] #Obtenemos los datos en ºK
lat = slp.variables[‘lat’][:] #Obtenemos los datos en º
lon = slp.variables[‘lon’][:] #Obtenemos los datos en º
[/sourcecode]
Los datos se han guardado en las variables slpdata, tsfcdata, lat y lon como numpy arrays. Para que la presión esté en hPa o mb (hectopascales o milibares), la temperatura en ºC y las longitudes vayan en una escala de -180º a 180º hacemos lo siguiente:
[sourcecode language=”python”]
slpdata = slpdata * 0.01
tsfcdata = tsfcdata – 273.15
lon[lon > 180] = lon – 360.
[/sourcecode]
Ahora creamos una instancia a Basemap:
[sourcecode language=”python”]
m = bm.Basemap(llcrnrlon = -20,
llcrnrlat = 25,
urcrnrlon = 60,
urcrnrlat = 80,
projection = ‘mill’)
[/sourcecode]
En el anterior código hemos puesto lo siguiente:
| llcrnrlon longitud de la esquina inferior izquierda del dominio del mapa seleccionado.
| llcrnrlat latitud de la esquina inferior izquierda del dominio del mapa seleccionado.
| urcrnrlon longitud de la esquina superior derecha del dominio del mapa seleccionado.
| urcrnrlat latitud de la esquina superior derecha del dominio del mapa seleccionado.
Para projection hemos usado mill que será para una proyección Miller cylindrical. Si queréis ver todas las proyecciones disponibles podéis escribir:
[sourcecode language=”python”]
print bm.supported_projections
[/sourcecode]
Y veréis el siguiente output:
[sourcecode language=”python”]
aeqd Azimuthal Equidistant
poly Polyconic
gnom Gnomonic
moll Mollweide
tmerc Transverse Mercator
nplaea North-Polar Lambert Azimuthal
gall Gall Stereographic Cylindrical
mill Miller Cylindrical
merc Mercator
stere Stereographic
npstere North-Polar Stereographic
hammer Hammer
geos Geostationary
nsper Near-Sided Perspective
vandg van der Grinten
laea Lambert Azimuthal Equal Area
mbtfpq McBryde-Thomas Flat-Polar Quartic
sinu Sinusoidal
spstere South-Polar Stereographic
lcc Lambert Conformal
npaeqd North-Polar Azimuthal Equidistant
eqdc Equidistant Conic
cyl Cylindrical Equidistant
omerc Oblique Mercator
aea Albers Equal Area
spaeqd South-Polar Azimuthal Equidistant
ortho Orthographic
cass Cassini-Soldner
splaea South-Polar Lambert Azimuthal
robin Robinson
[/sourcecode]
Buscamos los valores x e y que representarán a las latitudes y longitudes de la proyección del mapa seleccionada:
[sourcecode language=”python”]
## Encontramos los valores x,y para el grid de la proyección del mapa.
lon, lat = np.meshgrid(lon, lat)
x, y = m(lon, lat)
[/sourcecode]
Vamos a representar el campo de presiones en el mapa:
[sourcecode language=”python”]
fig=plt.figure(figsize=(8,6))
ax = fig.add_axes([0.05,0.05,0.9,0.85])
cs = m.contour(x,y,slpdata[0,:,:],np.arange(900,1100.,5.),colors=’y’,linewidths=1.25)
m.drawparallels(np.arange(0,360,10),labels=[1,1,0,0])
m.drawmeridians(np.arange(-180,180,10),labels=[0,0,0,1])
m.bluemarble()
plt.show()
[/sourcecode]
Este último código mostrará la siguiente figura:
En el anterior mapa hemos mostrado las líneas de nivel de la presión, hemos dibujado meridianos y paralelos y lo hemos representado con una proyección cilíndrica de Miller usando como fondo los datos Blue Marble de la NASA.
Ahora vamos a introducir también los datos de temperatura usando contourf (la f viene de fill, relleno y son contornos rellenados). Metemos lo siguiente en nuestro script:
[sourcecode language=”python”]
# create figure.
fig=plt.figure(figsize=(8,6))
ax = fig.add_axes([0.05,0.05,0.9,0.85])
cs = m.contour(x,y,slpdata[0,:,:],np.arange(900,1100.,5.),colors=’k’,linewidths=1.)
csf = m.contourf(x,y,tsfcdata[0,:,:],np.arange(-50,50.,2.))
m.drawcoastlines(linewidth=1.25, color=’grey’)
m.drawparallels(np.arange(0,360,10),labels=[1,1,0,0])
m.drawmeridians(np.arange(-180,180,10),labels=[0,0,0,1])
plt.show()
[/sourcecode]
Y obtenemos el siguiente resultado:
Donde hemos representado la presión al nivel del mar (isolíneas en color negro), las temperaturas en superficie (contornos de color rellenos), las líneas de costa (líneas continuas grises), paralelos y meridianos.
El script final quedaría algo así:
[sourcecode language=”python”]
## Importamos las librerías que vamos a usar
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits import basemap as bm
## Abrimos los ficheros de datos
slp = nc.Dataset(‘X83.34.8.250.104.4.18.19.nc’) #slp por ‘sea level pressure’
tsfc = nc.Dataset(‘X83.34.8.250.104.4.15.31.nc’) #tsfc ‘por temperature at surface’slp.variables
print slp.variables
print tsfc.variables
slpdata = slp.variables[‘slp’][:] #Obtenemos los datos en Pa
tsfcdata = tsfc.variables[‘air’][:] #Obtenemos los datos en ºK
lat = slp.variables[‘lat’][:] #Obtenemos los datos en º
lon = slp.variables[‘lon’][:] #Obtenemos los datos en º
slpdata = slpdata * 0.01
tsfcdata = tsfcdata – 273.15
lon[lon > 180] = lon – 360.
## Creamos una instancia a Basemap
m = bm.Basemap(llcrnrlon = -20,
llcrnrlat = 25,
urcrnrlon = 60,
urcrnrlat = 80,
projection = ‘mill’)
print bm.supported_projections
## Encontramos los valores x,y para el grid de la proyección del mapa.
lon, lat = np.meshgrid(lon, lat)
x, y = m(lon, lat)
# Creamos la figura con P sobre fondo Blue Marble.
fig=plt.figure(figsize=(8,6))
ax = fig.add_axes([0.05,0.05,0.9,0.85])
cs = m.contour(x,y,slpdata[0,:,:],np.arange(900,1100.,5.),colors=’y’,linewidths=1.25)
m.drawparallels(np.arange(0,360,10),labels=[1,1,0,0])
m.drawmeridians(np.arange(-180,180,10),labels=[0,0,0,1])
m.bluemarble()
plt.show()
# Creamos la figura con P y T.
fig=plt.figure(figsize=(8,6))
ax = fig.add_axes([0.05,0.05,0.9,0.85])
cs = m.contour(x,y,slpdata[0,:,:],np.arange(900,1100.,5.),colors=’k’,linewidths=1.)
csf = m.contourf(x,y,tsfcdata[0,:,:],np.arange(-50,50.,2.))
m.drawcoastlines(linewidth=1.25, color=’grey’)
m.drawparallels(np.arange(0,360,10),labels=[1,1,0,0])
m.drawmeridians(np.arange(-180,180,10),labels=[0,0,0,1])
plt.show()
[/sourcecode]
Y eso es todo por hoy. En algún momento, siempre que el tiempo lo permita, veremos más en profundidad Basemap y Netcdf4-python.
Saludos.
P.D.: Lo de siempre, si encontráis errores, queréis criticar (constructivamente) mis bajas dotes como programador o queréis aportar alguna cosa usad los comentarios.
Bueno, los resultados me parecen alucinantes. ¡Muchas gracias por la entrada Kiko! Esto ha sido un ejemplo de compenetración sin precedentes 😛
Pingback: Regiones de estabilidad de métodos numéricos en Python « Pybonacci
Esta espectacular el tuto… a ver si me puesde ayudar… tengo varios archivos extensión nc con lat lon y concentración de ozono cada archivo corresponde a un año y dentro de el están datos para 365 días. Necesito extraer los datos y hacer una serie de tiempo para una lat y lon determinada… HELPPPP
Los datos son
ftp://outgoing:outgoing@ftp.bodekerscientific.com/CombinedTCOV2.8/Daily/LongPatched/NetCDF/
Hola.
Localiza a que valores del grid x, y, corresponde tu posición lon, lat y extrae solo los 365 valores de t. Si quieres la serie de todos los años haz un bucle sobre todos los ficheros ordenados por año.
Pingback: ¿Por qué usar netCDF? « Pybonacci
Gracias por tu ayuda, pero me gustaria que me ayudarás con un archivo nc que cuenta con datoa desde 1900 hasta 2011, pero necesitmo datos desde 1981 hasta 2010, como haría para seleccionar esos datos.
Pingback: Basemap y Google Geocode para representar puntos sobre un mapa | Pybonacci
excelentes gente!! saben hacia falta que unos genios se pusieran a inventar un blog sobre la importancia que tiene la programacion cientifica y de hecho para estudiantes como yo (meteorologia ) quede fascinado con este articulo.. saludos y sigan asi
Gracias por compartir este tutorial, trabajo con WRF y pensaba en sacar los resultados via ncl. Definitivamente con Python, todo es más sencillo! Un saludo
Yo también trabajo con wrf y ncl/grads/otros no son una gran solución. Mejor depender de un lenguaje de programación completo y maravilloso como Python 😛
Hola Kiko. Los enlaces para descargar los ficheros .nc no están funcionando. Me enviasn a > “Oops!
Sorry, there is a problem with the server configuration.”
Yo quisiera hacer algo similar con los datos de viento sobre el océano, generados por el ASCAT… sabes si hay algo sobre eso_? o que recomendación me darías para hacerlo yo mismo?
Gracias Kiko~
Parece que han cambiado los scripts de descarga, luego lo miro en casa.
Para descargar datos ASCAT puedes echarle un ojo aquí: http://podaac.jpl.nasa.gov/datasetlist?search=ascat
Los tienes disponibles en netcdf comprimidos. Además, en esa misma página tienes también viento a partir de otros sensores como el QuikScat com un periodo de datos mucho más largo que el del ASCAT.
Todo se puede descargar desde ftp, openDAP, THREDDS,…, más fácil imposible,
Si tienes algún problema con los datos ASCAT puedes preguntar en http://es.stackoverflow.com y nos avisas por twitter o correo y te intentamos responder por ahí.
Saludos.
Muy amable, gracias!
Hola Kiko! Excelente tutorial, me ha ayudado mucho. Tengo una duda, si quiero graficar otro dia del año y no 01/01/2012. Cómo lo cambio en el código que has publicado?
Gracias!!!
Hola, Yania.
Tendrías que ir a la página de reanálisis del NCEP y descargar los días, variables, área y niveles de altura que quieras graficar.
http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html
Saludos.
P.D.: Si tienes problemas puedes consultar en http://es.stackoverflow.com/, que es un medio más apropiado para ello. Alguno de nosotros te contestará por ahí.