Grib-säädata graafiseksi esitykseksi kartalle

Ilmatieteen laitos julkaisee avoimena datana GRIB-muodossa joitakin sääparametreja. Grib on binäärimuodossa olevaa meteorologista dataa, joka koostuu selittävistä otsaketiedoista sekä itse säädatasta. Gribejä on saatavissa vanhemmassa grib1- ja uudemmassa grib2-muodossa. Jälkimmäinen on pakattavissa paremmin.

Koska tiedosto on binäärimuodossa, tarvitaan jokin ohjelma, joka osaa sitä lukea. Maksullisia ohjelmia on monia, mutta maksuttomia vähemmän. Olen aiemmin kuvannut Grads-ohjelmaa, jolla voi ilmaiseksi katsella grib-tiedoston sisältöä graafisesti. Grads on kuitenkin melko rajoittunut, tai ainakin osaamiseni sen kanssa on rajoittunut, joten tässä artikkelissa käsitellään python-kirjastoja, joilla saadaan grib-data esitettyä graafisesti kartalla.

Pygrib on python-ohjelmiontikielen moduuli, joka osaa lukea grib-tiedostoja (sekä 1 että 2). Pygrib hyödyntää Eurooppakeskuksen (ECMWF) julkaisemaa c-kirjastoa GRIB APIa. Seuraavassa esitellään asennus Ubuntussa, vaikka samat asiat onnistunevat soveltaen myös muissa ympäristöissä.

Pygribin asennusta varten tarvitaan

  1. python 2.4 tai uudempi, itsellä Python 2.7.4 (sudo apt-get install python2.7)
  2. numpy python-laajennus N-ulotteinille taulukko-olioille, 1.2.1 tai uudempi (sudo apt-get install python-numpy)
  3. matplotlib kirjasto graafisten esitysten piirtoon (sudo apt-get install python-matplotlib)
  4. basemap toolkit (sudo apt-get install python-mpltoolkits.basemap)
  5. Jasper jpeg-tiedostoille
  6. libpng png-tiedostoille (sudo apt-get install libpng-dev)
  7. GRIB API c-kirjasto grib-tiedostojen aukaisemiseen ja tallentamiseen (Ubuntulle löytyy oma .deb-paketti)’
  8. Pygrib

    • export GRIBAPI_DIR=polku_grip_apiim
    • export JASPER_DIR=polku_jasperiin
    • python setup.py install –user
    • python test.py (jos kaikki meni hyvin, virheitä ei pitäisi ilmetä)

Kun pygrib on asennettu, voidaan sitä ja muita kirjastoja apuna käyttäen muodostaa esimerkiksi seuraavanlainen ohjelma, joka lukee grib-tiedoston, piirtää kartan sekä sen päälle täytetyt lämpötilakontuurit. Ohjelman ymmärtämiseksi on syytä opetella ainakin basemapin ja matplotlibin toimintaa.

import pygrib, time
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
 
grbs = pygrib.open('hirlam.grb')
btemps = [grb for grb in grbs if grb['name']=='2 metre temperature']
grb = btemps[0]
grbs.close()
 
m = Basemap(llcrnrlat=59.24,llcrnrlon=18.8,urcrnrlat=70.45,urcrnrlon=35.03,
            resolution='l',projection='tmerc',lat_0=65.091,lon_0=27.24)
 
fig = plt.figure(figsize=(8,7))
 
for i, cp in enumerate(m.coastpolygons):
     if m.coastpolygontypes[i] < 2:
         m.plot(cp[0],cp[1],'k-', zorder=4, color='black', linewidth=1, antialiased=1)
 
m.fillcontinents(color='coral', lake_color='coral')
m.drawmapboundary(fill_color='aqua')
 
ny = grb['values'].shape[0]
nx = grb['values'].shape[1]
lons, lats = m.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid.
x, y = m(lons, lats) # compute map proj coordinates.
# draw filled contours.
clevs = [255,260,265,270,275,280,285,290,295]
 
cs = m.contourf(x, y, grb['values'], clevs, alpha=1, antialiased=True, cmap=plt.cm.rainbow)
 
for col in cs.collections:
     col.set_zorder(2.5) 
m.drawcountries(linewidth=1, linestyle='solid', color='black', antialiased=1, ax=None, zorder=3)
 
plt.show()

Tällä saadaan seuraavanlainen kuva lämpötilakontuureista (Kelvineinä gribissä):

finland_temperature_contour

Tässä tosin data ja kartta eivät osu aivan kohdilleen eri projektioiden vuoksi. Datassa käytetään EPSG:4326-projektiota ja kartalle on määritetty transverse mercator. Alapuolella on kartan ja datan projektiot samat, mutta tällöin Suomi näyttää lättänältä:

finland_cylindrical_equidistant

Eli datan projektiota pitäisi muokata…

Päivitys 15.10.2013

Huomasinpa tuossa, että basemap onkin sen verran viisas, että osaa hoitaa datan projisoinnin ihan itsestäänkin, kunhan vain antaa sille kaikki datan koordinaatit. Eli edellisessä koodissa pitää korvata rivit

ny = grb['values'].shape[0]
nx = grb['values'].shape[1]
lons, lats = m.makegrid(nx, ny)

seuraavalla

lats, lons = grb.latlons()

ja nyt tuleekin jo ihan suht järkevän näköistä outputia:

2_metre_temperature_002

Julkaistu kohteessa Uncategorized

Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *

*