ecmwf / cfgrib

A Python interface to map GRIB files to the NetCDF Common Data Model following the CF Convention using ecCodes
Apache License 2.0
400 stars 77 forks source link

Special behavior with wind and temperature ? #249

Closed abakleriche closed 3 years ago

abakleriche commented 3 years ago

I have GRIB files with parameters at 25 levels of height. If I load the file of the geopotential, I get a dataset with a 4D variable Z But if I load the U component of wind, I only get a 3D variable u10 (wind at 10 m). U

The other levels are considered as another variable, which I can retrieve in another dataset if I use open_datasets instead of open_dataset. I observe the same behavior with the temperature at 2 m. The problem is that if I use xarray.open_mfdataset to open U file and V file in order to get them in a unique dataset, I only get u10 and v10. Is there a means to cancel this special behavior for U, V and T ?

iainrussell commented 3 years ago

Hello @abakleriche - would you be able to attach these GRIB files so that we can have a closer look please? I'm guessing that cfgrib is getting confused between surface and height-level fields.

Thanks, Iain

abakleriche commented 3 years ago

Hi, Thank you for answering. Unfortunately these files are quite big (more than 40 Mb each after zipping), so I think I can't upload them in Github. When I try I get a message "is not included to the list". Is there another way to upload them ?

iainrussell commented 3 years ago

Hi @abakleriche , Could you upload to Dropbox or similar cloud storage?

Alternatively, perhaps you have tools that could help you to reduce the size of the files, for example by reducing the bits per value, e.g. see the last example here: https://confluence.ecmwf.int/display/ECC/grib_set or this: https://confluence.ecmwf.int/display/UDOC/How+can+I+remove+the+data+values+and+bitmap+-+ecCodes+GRIB+FAQ

Cheers, Iain

abakleriche commented 3 years ago

Hi Iain,

I uploaded the files to this link : https://www.transfernow.net/dl/20210920qeVqCg1K (valid for 7 days) Best regards,

Roger

iainrussell commented 3 years ago

Hi Roger,

Many thanks for uploading the files. I can see the cause of the problem now, though I'm not quite sure what the solution will be. The problem comes from the fact that the ecCodes library has a special understanding of some of these fields. In particular 10m wind components are identified specifically as '10U' and '10V', whereas the other height levels of wind are just called 'U' and 'V'. This can be seen with grib_ls command on the command line, or using Metview's Grib Examiner application or Metview's Python ls() command:

import metview
u = metview.read('TEMPETE_HAUTEUR_2017030100_U.grib')
u.ls()
---
        centre shortName    typeOfLevel      level  dataDate  dataTime stepRange dataType number   gridType 
Message                                                                                                     
0        lfpw      10u    heightAboveGround    10   20170301      0         0      None    None   regular_ll
1        lfpw        u    heightAboveGround    20   20170301      0         0      None    None   regular_ll
2        lfpw        u    heightAboveGround    35   20170301      0         0      None    None   regular_ll
3        lfpw        u    heightAboveGround    50   20170301      0         0      None    None   regular_ll

This causes cfgrib some problems because they are seen as different parameters on different sets of levels (10U is on level 10, whereas U is on levels 20 to 3000).

One thing you can do is to create two xarray datasets by applying a filter, e.g.

dsu = xr.open_dataset('TEMPETE_HAUTEUR_2017030100_U.grib',
                      engine='cfgrib',
                      backend_kwargs={'filter_by_keys': {'shortName': 'u'}}
                      )
dsu10 = xr.open_dataset('TEMPETE_HAUTEUR_2017030100_U.grib',
                      engine='cfgrib',
                      backend_kwargs={'filter_by_keys': {'shortName': '10u'}}
                      )

Now, it might be possible through some xarray trickery to combine these into a single dataset. I'm not the expert in doing that, but I suspect it should be possible.

You'll see the same issue with V vs 10V and with T vs 2T. I hope this helps you to proceed though!

Best regards, Iain

abakleriche commented 3 years ago

I Iain,

Thank you for your answer.

Best regards,

Roger