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:
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)
source to share
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.).
source to share
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!
source to share