Search code examples
pythonmatplotlibmatplotlib-basemapgribgfs

White area on matplotlib plot with pygrib data between 359.5 and 360 degrees


What I try is to plot output of the gfs weather model with matplotlib using pygrib to save the data, which is saved in grib files. Nearly everything works fine, the output looks like this:

enter image description here

It appears that the program isn't closing the gap between 359.5 and 360 degress by using the data of 0 degress. If the data would be in a regular list or something I would use the data of 0° and save it for 360° too by appending the list. I've seen people having the same problem with non-pygrib data. If you know how to change the pygrib data (regular operations don't work on pygrib data unfortunately) or how to make matplotlib close the gap, you would really help me out of this problem. Maybe the function "addcyclic" could help, but I don't know how.

EDIT: I solved the problem, see my answer.

So here is the code producing the problem:

#!/usr/bin/python3

import os, sys, datetime, string
from abc import ABCMeta, abstractmethod

import numpy as np
import numpy.ma as ma
from scipy.ndimage.filters import minimum_filter, maximum_filter
import pygrib
from netCDF4 import Dataset
from pylab import *
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid

import laplaceFilter
import mpl_util


class Plot(Basemap):
    def __init__(self, basemapParams):
        super().__init__(**basemapParams)
        self.layers = []

    def addLayer(self, layer):
        self.layers.append(layer)

    def plot(self, data):
        for layer in self.layers:
            layer.plot(self, data)

        plt.title('Plot')
        plt.show()


class Layer(metaclass=ABCMeta):
    def __init__(self):
        pass

    @abstractmethod
    def plot(self, plot, data):
        return NotImplemented


class BackgroundLayer(Layer):
    def __init__(self, bgtype, coords):
        #possible bgtype values: borders, topo, both
        self.bgtype = bgtype
        self.lonStart = coords[0]
        self.lonEnd   = coords[1]
        self.latStart = coords[2]
        self.latEnd   = coords[3]

    def plot(self, plot, data):
        [...]


    def findSubsetIndices(self,min_lat,max_lat,min_lon,max_lon,lats,lons):
        [...]


class LegendLayer(Layer):
    def __init__(self):
        pass


class GribDataLayer(Layer, metaclass=ABCMeta):
    def __init__(self, varname, level, clevs, cmap, factor):
        self.varname = varname
        self.level = level
        self.clevs = clevs
        self.cmap = cmap
        self.factor = factor

    def plot(self, plot, data):
        #depending on the height we want to use, we have to change the index
        indexes = {1000:0, 2000:1, 3000:2, 5000:3, 7000:4, 10000:5, 15000:6, 20000:7, 25000:8, 30000:9,
                   35000:10, 40000:11, 45000:12, 50000:13, 55000:14, 60000:15, 65000:16, 70000:17,
                   75000:18, 80000:19, 85000:20, 90000:21, 92500:22, 95000:23, 97500:24, 100000:25, 0:0}

        selecteddata = data.select(name = self.varname)[indexes[self.level]]
        lats, lons = selecteddata.latlons()



        layerdata = selecteddata.values*self.factor

        x, y = plot(lons, lats) # compute map proj coordinates.
        self.fillLayer(plot, x, y, layerdata, self.clevs, self.cmap)

    @abstractmethod
    def fillLayer(self, plot, x, y, layerdata, clevs, cmap):
        return NotImplemented


class ContourLayer(GribDataLayer):
    def __init__(self, varname, level, clevs, cmap, factor, linewidth=1.5, fontsize=15,
                 fmt="%3.1f", inline=0,labelcolor = 'k'):
        self.linewidth = linewidth
        self.fontsize = fontsize
        self.fmt = fmt
        self.inline = inline
        self.labelcolor = labelcolor
        super().__init__(varname, level, clevs, cmap, factor)

    def fillLayer(self, plot, x, y, layerdata, clevs, cmap):
        # contour data over the map.
        cs = plot.contour(x,y,layerdata,clevs,colors = cmap,linewidths = self.linewidth)
        plt.clabel(cs, clevs, fontsize = self.fontsize, fmt = self.fmt, 
                   inline = self.inline, colors = self.labelcolor)
        if self.varname == "Pressure reduced to MSL":
            self.plotHighsLows(plot,layerdata,x,y)

    def plotHighsLows(self,plot,layerdata,x,y):
        [...]


class ContourFilledLayer(GribDataLayer):
    def __init__(self, varname, level, clevs, cmap, factor, extend="both"):
        self.extend = extend
        super().__init__(varname, level, clevs, cmap, factor)

    def fillLayer(self, plot, x, y, layerdata, clevs, cmap):
        # contourfilled data over the map.
        cs = plot.contourf(x,y,layerdata,levels=clevs,cmap=cmap,extend=self.extend)
        #cbar = plot.colorbar.ColorbarBase(cs)


[...]


ger_coords = [4.,17.,46.,56.]
eu_coords  = [-25.,57.,22.,70.]


### Choose Data
data = pygrib.open('gfs.t12z.mastergrb2f03')


### 500hPa Europe
coords = eu_coords
plot1 = Plot({"projection":"lcc","resolution":"h","rsphere":(6378137.00,6356752.3142), "area_thresh": 1000.,
             "llcrnrlon":coords[0],"llcrnrlat":coords[2],"urcrnrlon":coords[1],"urcrnrlat":coords[3],
             "lon_0":(coords[0]+coords[1])/2.,"lat_0":(coords[2]+coords[3])/2.})


clevs = range(480,600,4)
cmap = plt.cm.nipy_spectral
factor = .1
extend = "both"
level = 50000
layer1 = ContourFilledLayer('Geopotential Height', level, clevs, cmap, factor, extend)


clevs = [480.,552.,600.]
linewidth = 2.
fontsize = 14
fmt = "%d"
inline = 0
labelcolor = 'k'
layer2 = ContourLayer('Geopotential Height', level, clevs, 'k', factor, linewidth, fontsize, fmt, inline, labelcolor)

level = 0
clevs = range(800,1100,5)
factor = .01
linewidth = 1.5
inline = 0
labelcolor = 'k'
layer3 = ContourLayer('Pressure reduced to MSL', level, clevs, 'w', factor, linewidth, fontsize, fmt, inline, labelcolor)


plot1.addLayer(BackgroundLayer('borders', coords))
plot1.addLayer(layer1)
plot1.addLayer(layer2)
plot1.addLayer(layer3)
plot1.plot(data)

Solution

  • I solved it myself 2 months later:

    Matplotlib doesn't fill the area if your longitude range is from 0 to 359.75 because it ends there from matplotlibs point of view. I solved it by dividing up the data and then stacking it.

    selecteddata_all = data.select(name = "Temperature")[0]
    
    selecteddata1, lats1, lons1 = selecteddata_all.data(lat1=20,lat2=60,lon1=335,lon2=360)
    selecteddata2, lats2, lons2 = selecteddata_all.data(lat1=20,lat2=60,lon1=0,lon2=30)
    
    lons         = np.hstack((lons1,lons2))
    lats         = np.hstack((lats1,lats2))
    selecteddata = np.hstack((selecteddata1,selecteddata2))
    

    No white area left of 0° anymore.

    I don't know whether there is a fix if you wanna plot a whole hemisphere (0 to 359.75 deg).