matplotlib / basemap

Plot on map projections (with coastlines and political boundaries) using matplotlib
MIT License
774 stars 390 forks source link

addcyclic problem with grib2 data #287

Open fox91 opened 8 years ago

fox91 commented 8 years ago

I have some problems to use add cyclic with data from grib2 file.

This is my code:

# -*- coding: utf-8 -*-

import pygrib
import numpy as np

from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure

from mpl_toolkits.basemap import Basemap, addcyclic

fig = Figure((4, 3))
canvas = FigureCanvas(fig)
ax = fig.add_axes([0.05, 0.05, 0.9, 0.9])

# File download from http://para.nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/para/gfs.2016041918/gfs.t18z.pgrb2.0p25.f000
grbindx = pygrib.index('gfs.t18z.pgrb2.0p25.f000', 'shortName', 'typeOfLevel', 'level')
grb = grbindx.select(shortName='2t', typeOfLevel='heightAboveGround', level=2)[0]
lats, lons  = grb.latlons()
temp = grb.values

temp, lons = addcyclic(temp, lons)

bm = Basemap(projection='mill', resolution='l', ax=ax,
    llcrnrlat=32, urcrnrlat=70, llcrnrlon=-25,
    urcrnrlon=40.5, lat_ts=20)

x, y = bm(lons, lats)
bm.contourf(x, y, temp)

bm.drawcoastlines()

canvas.print_figure('gfs_temperature_2m.png', dpi=100)

With python 2.7.11 and basemap 1.0.7 this is the result:

Traceback (most recent call last):
  File "addcyclic.py", line 21, in <module>
    temp, lons = addcyclic(temp, lons)
  File "/usr/local/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py", line 5062, in addcyclic
    lonsout[0:nlons] = lonsin[:]
ValueError: could not broadcast input array from shape (721,1440) into shape (1440)

With python 2.7.11 and basemap 1.0.8 (from latest github commit) this is the result:

Traceback (most recent call last):
  File "addcyclic.py", line 27, in <module>
    x, y = bm(lons, lats)
  File "/usr/local/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.py", line 1175, in __call__
    xout,yout = self.projtran(x,y,inverse=inverse)
  File "/usr/local/lib/python2.7/site-packages/mpl_toolkits/basemap/proj.py", line 286, in __call__
    outx,outy = self._proj4(x, y, inverse=inverse)
  File "/usr/local/lib/python2.7/site-packages/pyproj/__init__.py", line 399, in __call__
    _proj.Proj._fwd(self, inx, iny, radians=radians, errcheck=errcheck)
  File "_proj.pyx", line 128, in _proj.Proj._fwd (_proj.c:1678)
RuntimeError: Buffer lengths not the same

Where is the problem? Some bug in basemap? Or in my code?

Without addcycle this is the result: gfs_temperature_2m

Thanks Andrew

guziy commented 8 years ago

Maybe the shapes of the data and lons arrays should be the same, try using np.meshgrid:

lons2d, lats2d = np.meshgrid(lons, lats)
# use lons2d for addcyclic... and plotting

Cheers

fox91 commented 8 years ago

Hi @guziy, your line of code give me this error:

Traceback (most recent call last):
  File "addcyclic.py", line 21, in <module>
    lons2d, lats2d = np.meshgrid(lons, lats)
  File "/usr/local/lib/python2.7/site-packages/numpy/lib/function_base.py", line 4104, in meshgrid
    mult_fact = np.ones(shape, dtype=int)
  File "/usr/local/lib/python2.7/site-packages/numpy/core/numeric.py", line 190, in ones
    a = empty(shape, dtype, order)
MemoryError
guziy commented 8 years ago

At this point, I would suggest uploading your file somewhere, just to check if meshgrid solves your initial problem.

Cheers

fox91 commented 8 years ago

@guziy are you referring to GRIB2 data file? You can download it from http://para.nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/para/gfs.YYYYMMDD18/gfs.t18z.pgrb2.0p25.f000 using yesterday date instead of YYYYMMDD. Instead my python script is here.

I have this python packages installer:

basemap (1.0.7)
cycler (0.10.0)
Cython (0.24)
h5py (2.5.0)
matplotlib (1.5.1)
netCDF4 (1.2.3.1)
numpy (1.11.0)
OWSLib (0.11.0)
Pillow (3.2.0)
pip (8.1.1)
pygrib (2.0.1)
pyparsing (2.1.1)
pyproj (1.9.5.1)
pyshp (1.2.3)
python-dateutil (2.5.2)
pytz (2016.3)
requests (2.9.1)
scipy (0.17.0)
setuptools (20.6.7)
six (1.10.0)
virtualenv (15.0.1)
wheel (0.29.0)

Thanks

guziy commented 8 years ago

meshgrid is failing because your lons and lats are already 2d.

fox91 commented 8 years ago

@guziy yes, in fact, in my example I had not included. I upgraded my gist.

But my question is: how do I properly use addcyclic with this data? The problem does not occur only with this file, but with all GRIB2 files. On the contrary, by reading data from NetCDF I have no problems.

Thanks

guziy commented 8 years ago

If you need the specified region you would have to shift the part of data and longitudes (those that are greater than 180) and make sure that longitudes increase from left to right.

I actually was trying to use shiftgrid to shift your data and longitudes, but that did not help... I'll let you know if I come up with something better to solve this.

Cheers

fox91 commented 8 years ago

@guziy thanks for your help. In this example I centered the map on Europe to show the problem. I also need all world map, but from -180 to 180. Instead data start from 0 to 360. addcyclic should not be used for this?

Cheers

guziy commented 8 years ago

I think the shiftgrid function is used for what you need. But unfortunately, it assumes 1d coordinates...

The most general way of doing it would be to convert your coordinates to cartesian with the center at the Earth center, and then use nearest neighbor interpolation to the target grid.

You could check this post for interpolation: http://earthpy.org/interpolation_between_grids_with_ckdtree.html

There is another post, which uses pyresample, worth checking out.

Cheers, sorry for the late response, soo many things going on))