Python でお絵かき

目的

  • GPhys, GGraph チュートリアルでやっていることを Python でやろうとするとどうなるかやってみた

準備すること

Python の環境

  • anaconda でインストール (Python 3.6)
  • NetCDF4 を読むための環境を整える(以下、わたしの作業ログ)
  $ conda install -c anaconda netcdf4=1.2.4    # NetCDF4 モジュール
  $ conda install -c anaconda hdf4             # libmfhdf.so.0 がないといわれたので
  $ conda install -c anaconda basemap=1.0.7    # Basemap モジュール
In [5]:
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
data = Dataset("air.2012-01.nc", "r")
print(data)
<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_CLASSIC data model, file format NETCDF3):
    Conventions: COARDS
    title: mean daily NMC reanalysis (2012)
    history: created 2011/12 by Hoop (netCDF2.3)
2012-02-28 18:40:29 JST horinout> extractJan.rb wrote air
    description: Data is from NMC initialized reanalysis
(4x/day).  It consists of most variables interpolated to
pressure surfaces from model (sigma) surfaces.
    platform: Model
    references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html
    dimensions(sizes): lon(144), lat(73), level(17), time(31)
    variables(dimensions): float32 lon(lon), float32 lat(lat), float32 level(level), float64 time(time), int16 air(time,level,lat,lon)
    groups: 

In [6]:
data.variables.keys()
Out[6]:
odict_keys(['lon', 'lat', 'level', 'time', 'air'])
In [7]:
data.variables['lon']
Out[7]:
<class 'netCDF4._netCDF4.Variable'>
float32 lon(lon)
    units: degrees_east
    long_name: Longitude
    actual_range: [   0.   357.5]
    standard_name: longitude
    axis: X
unlimited dimensions: 
current shape = (144,)
filling off
In [8]:
lon = data.variables['lon'][:]
print(lon)
[   0.     2.5    5.     7.5   10.    12.5   15.    17.5   20.    22.5
   25.    27.5   30.    32.5   35.    37.5   40.    42.5   45.    47.5
   50.    52.5   55.    57.5   60.    62.5   65.    67.5   70.    72.5
   75.    77.5   80.    82.5   85.    87.5   90.    92.5   95.    97.5
  100.   102.5  105.   107.5  110.   112.5  115.   117.5  120.   122.5
  125.   127.5  130.   132.5  135.   137.5  140.   142.5  145.   147.5
  150.   152.5  155.   157.5  160.   162.5  165.   167.5  170.   172.5
  175.   177.5  180.   182.5  185.   187.5  190.   192.5  195.   197.5
  200.   202.5  205.   207.5  210.   212.5  215.   217.5  220.   222.5
  225.   227.5  230.   232.5  235.   237.5  240.   242.5  245.   247.5
  250.   252.5  255.   257.5  260.   262.5  265.   267.5  270.   272.5
  275.   277.5  280.   282.5  285.   287.5  290.   292.5  295.   297.5
  300.   302.5  305.   307.5  310.   312.5  315.   317.5  320.   322.5
  325.   327.5  330.   332.5  335.   337.5  340.   342.5  345.   347.5
  350.   352.5  355.   357.5]
In [9]:
lat = data.variables['lat'][:]
air = data.variables['air'][:]
print(air.shape)
(31, 17, 73, 144)
In [10]:
from mpl_toolkits.basemap import Basemap

plt.figure(figsize=(8,6))
bm = Basemap(projection="cyl", # Equidistant Cylindrical projection
            llcrnrlon=0, llcrnrlat=-70,
            urcrnrlon=360, urcrnrlat=70)
# labels = [left,right,top,bottom]
bm.drawcoastlines(linewidth=1.0) # 海岸線
bm.drawmeridians(np.arange(0, 360, 30), labels=[0,0,0,1]) # 経線
bm.drawparallels(np.arange(-90, 90, 30), labels=[1,0,0,0]) # 緯線
Out[10]:
{-60: ([<matplotlib.lines.Line2D at 0x7f157b237860>],
  [<matplotlib.text.Text at 0x7f157b242630>]),
 -30: ([<matplotlib.lines.Line2D at 0x7f157b5145f8>],
  [<matplotlib.text.Text at 0x7f157b242b00>]),
 0: ([<matplotlib.lines.Line2D at 0x7f157b514dd8>],
  [<matplotlib.text.Text at 0x7f157b257080>]),
 30: ([<matplotlib.lines.Line2D at 0x7f157b4275f8>],
  [<matplotlib.text.Text at 0x7f157b2575c0>]),
 60: ([<matplotlib.lines.Line2D at 0x7f157b427dd8>],
  [<matplotlib.text.Text at 0x7f157b257b00>])}
In [11]:
X, Y = np.meshgrid(lon, lat)
x, y = bm(X, Y)
cs = bm.pcolor(x, y, air[0,0,:,:], cmap='jet')
cbar = bm.colorbar(cs, location='right', pad='10%')
plt.show()
/home/matsuba/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3342: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.
  b = ax.ishold()
/home/matsuba/anaconda3/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py:3381: MatplotlibDeprecationWarning: axes.hold is deprecated.
    See the API Changes document (http://matplotlib.org/api/api_changes.html)
    for more details.
  ax.hold(b)
In [12]:
!gplist air.2012-01.nc
air.2012-01.nc:
  lon	[lon=144]	'Longitude'	(degrees_east)
  lat	[lat=73]	'Latitude'	(degrees_north)
  level	[level=17]	'Level'	(millibar)
  time	[time=31]	'Time'	(hours since 1-1-1 00:00:0.0)
  air	[lon=144,lat=73,level=17,time=31]	'mean Daily Air temperature'	(degK)

gpview で見てみるとどうなるか

In [13]:
!gpview air.2012-01.nc@air --nocont --itr 10 --map coast_world --clrmap 15
 *** MESSAGE (SWDOPN) ***  GRPH1 : STARTED / IWS =  1.                         
  Reading air.2012-01.nc@air
 *** WARNING (SLZTTL) ***  SPACE FOR TITLE IS NOT ENOUGH.                      
 *** MESSAGE (-CNT.-) ***  PY = 0.0 IS ASSUMED.                                
 *** WARNING (SLZTTL) ***  SPACE FOR TITLE IS NOT ENOUGH.                      
 *** MESSAGE (-CNT.-) ***  PY = 0.0 IS ASSUMED.                                
 *** MESSAGE (SWPCLS) ***  GRPH1 : PAGE =   1 COMPLETED.                       
 *** MESSAGE (SWDCLS) ***  GRPH1 : TERMINATED.                                 

ウィンドウが出現するので d ボタンで画像ダンプを作成する。その結果が以下である。

比較しやすいように図のサイズが自在に変更できればよいが、勉強不足なので、いまはこれでご容赦。