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

What I am trying to do is to display the output of the gfs weather model with matplotlib using pygrib to save the data that is saved in grib files. Almost everything works fine, the output looks like this:

enter image description here

It looks like the program doesn't close the gap between 359.5 and 360 degress using 0 degress data. If the data is going to be in a regular list or something, I will use 0 ° data and keep it 360 ° by adding the list. I have seen people having the same problem with non pygrib data. If you know how to change pygrib data (normal operations don't work with pygrib data, unfortunately) or how to get matplotlib to close the gap, you can really help me with this problem. Maybe the "addcyclic" function can help, but I don't know how.

EDIT: I solved the problem, see my answer.

So, here is the code creating 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)

      

+3


source to share


2 answers


I solved this myself after 2 months:

Matplotlib does not fill the area if your longitude range is 0 to 359.75 because it ends up from matplotlibs perspective. I solved this by splitting the data and 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 remaining from 0 ° anymore.

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

+1


source


I worked it myself myself and the base module add-on function really works really well. In the basemap documentation, lay out the syntax and use it pretty well.

In terms of variables in your code, you can add a looping dot either before or after you multiply self.factor in your GribDataLayer class:

layerdata, lons = addcyclic(layerdata, lons)

      

You can also use np.append and write your own function to accomplish the same task. It will look something like this:



layerdata = np.append(layerdata,layerdata[...,0,None],axis=-1)

      

If your input is 2D, then the above syntax is equivalent to fetching all data in the first strip of longitude (ie layerdata [:, 0])

layerdata = np.append(layerdata,layerdata[:,0,None],axis=-1)

      

Hope this helps!

0


source







All Articles