Visualize grib2 on a map with python (matplotlib)

Last time explained how to build an environment to operate grib2 files from python. This time the read data. I will introduce the method of visualization with matplotlib + basemap. I won't explain Jupyter in particular, but I'm basically trying it on Jupyter (I think it's okay to use python commands as usual).

Installation

In Anaconda environment

apt-get -y install libgl-dev
conda install matplotlib
conda install basemap

You can install matplotlib and Basemap with. The place of apt-get changes depending on the environment, but if you do not have it, you will be told that libGL.so does not exist when importing matplotlib. In the case of docker, it is OK if you add RUN to the beginning of each line above to the previous Dockerfile.

Dockerfile

...Add the following...
RUN apt-get -y install libgl-dev
RUN conda install matplotlib
RUB conda install basemap

grib2 data read

Use the pygrib described last time to read the grib data. Use the same sample as last time.

wget -P somepath/ http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2017/01/02/Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin

From here on, we'll work within python. When working with docker, start container below and start python in it as before.

cd path_of_Dokcerfile
docker build -t pygrib_ubuntu .
docker run -it -v /somepath:/temp pygrib_ubuntu

import matplotlib
matplotlib.use('Agg') #If there is no x environment on the server, matplotlib will not work without it
import pygrib
import matplotlib.pyplot as plt
import numpy as np
import math
from mpl_toolkits.basemap import Basemap


grbs = pygrib.open('/temp/Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')

#Element specified by select method(forecastTime=0)An array of elements that match is returned.
# forecastTime=There are multiple 0 elements, but the beginning([0])Is the sea level rehabilitation pressure(mslp)Is included
grb = grbs.select(forecastTime=0)[0]

# lats,lons is a two-dimensional array containing latitude and longitude
lats, lons = grb.latlons()

#Data contained in(Sea level rehabilitation pressure)With a two-dimensional ndarray
mslp = grb.values

With this, the latitude / longitude and the data at that point (sea level rehabilitation pressure) required for visualization can be extracted from the grib.

Visualization

Next, I will draw isobars like a weather map.

#List the intervals at which isobars are drawn in levels
#The unit of MSLP is Pa, so if you draw a line every 2hPa, you will draw a line every 200.
levels = range(math.floor(mslp.min()/100)*100, math.ceil(mslp.max()/100)*100+1,200)

# lat,Convert lon to one dimension
#For some reason it cannot be drawn in 2D
flat_lats= np.ravel(lats)
flat_lons= np.ravel(lons)

fig = plt.figure()

#Create a map. 4 parameters specify the drawing range
m = Basemap(llcrnrlat=lats.min(),urcrnrlat=lats.max(), llcrnrlon=lons.min(),urcrnrlon=lons.max())

#Draw isobars
m.contour(flat_lons, flat_lats,mslp,latlon=True,tri=True,levels=levels)

#Draw the coastline
m.drawcoastlines() 

fig.savefig('/temp/mslp.png')

The picture is completed like this.

The following warnings will appear at the time of output, but for the time being, it will not affect it, so let's not worry about it.

/root/.pyenv/versions/anaconda3-4.1.1/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3505: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.

  b = ax.ishold()

/root/.pyenv/versions/anaconda3-4.1.1/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3572: MatplotlibDeprecationWarning: axes.hold is deprecated.

    See the API Changes document (http://matplotlib.org/api/api_changes.html)

    for more details.

  ax.hold(b)