Ejemplo de uso de Basemap y NetCDF4

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]

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.

15 pensamientos sobre “Ejemplo de uso de Basemap y NetCDF4”

  1. Pingback: Regiones de estabilidad de métodos numéricos en Python « Pybonacci

  2. 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.

  3. Pingback: ¿Por qué usar netCDF? « Pybonacci

  4. 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.

  5. Pingback: Basemap y Google Geocode para representar puntos sobre un mapa | Pybonacci

  6. 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

  7. 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

    1. 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 😛

  8. 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~

    1. 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.

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

− three = 3