海面水温¶
地図だけでなく値を描画してみよう。ここでは海面水温データHadISST1を使う。
download pageからHadISST1_SST_update.nc.gzをダウンロードする。拡張子.gz（gzip圧縮）が付いているが，圧縮されていない場合は，拡張子.gz削除する。圧縮されている場合は解凍する。Windowsにはgzipを扱う機能はないので，7-Zipなどを用いる。
NetCDF¶
nc（NetCDF）を読むにはnetCDF4 pythonを用いる。
In [1]:
import netCDF4 as nc4
dataset = nc4.Dataset('HadISST1_SST_update.nc')
netcdf_dataset(ファイル名)はDatasetを作成する。Datasetは座標軸や変数などをまとめたものである。Datasetから変数を読むにはvariablesメソッドを用いる。
In [2]:
sst = dataset.variables['sst']
1月のデータを描画してみる。
In [3]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
lat = dataset.variables['latitude'][:]
lon = dataset.variables['longitude'][:]
ax = plt.axes(projection=ccrs.PlateCarree())
plt.contourf(lon, lat, sst[0,], 60, transform=ccrs.PlateCarree())
plt.colorbar()
ax.coastlines()
plt.show()
おかしな図になってしまった。データ形式を確認すると，海氷に覆われているところに-1000が入っている。確認してみる。
In [4]:
sst[0,:,180]
Out[4]:
masked_array(data = [-1000.0 -1000.0 -1.7999999523162842 -1.7999999523162842
-1.7999999523162842 -1.7999999523162842 -1.7999999523162842
-1.7999999523162842 -1.7999999523162842 -1.3447265625 -0.15625
1.0224127769470215 1.7339876890182495 2.135985851287842 1.9301509857177734
1.3141674995422363 1.8001981973648071 2.6573405265808105 3.610403299331665
4.4060893058776855 5.157680988311768 5.751012802124023 6.087249755859375
6.243302345275879 6.439548492431641 6.810280799865723 7.357912063598633
7.942025184631348 8.291414260864258 8.346809387207031 8.331205368041992
8.400185585021973 8.414026260375977 8.284310340881348 8.080513000488281
7.885481357574463 7.802855491638184 -- -- 9.766984939575195
9.766984939575195 -- -- -- -- -- -- -- -- 14.950439453125
15.034799575805664 15.309090614318848 15.588996887207031
15.678974151611328 -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- 28.66602897644043 28.807947158813477
28.998897552490234 28.993816375732422 28.805492401123047 28.4868221282959
28.14432716369629 27.855770111083984 27.61625099182129 27.427505493164062
27.250234603881836 27.02727508544922 26.753498077392578 26.424318313598633
26.01759910583496 25.558128356933594 25.080873489379883 24.617101669311523
24.203371047973633 23.833818435668945 23.51828384399414 23.262577056884766
23.099180221557617 23.022619247436523 22.979427337646484 22.97173500061035
22.976743698120117 22.99199867248535 23.07520866394043 23.2158145904541
23.322978973388672 23.351940155029297 23.316436767578125
23.211238861083984 23.045005798339844 22.852628707885742
22.544879913330078 22.095966339111328 21.5771484375 20.952518463134766
20.155269622802734 19.183979034423828 18.14966583251953 17.092430114746094
16.00191307067871 14.901739120483398 13.715605735778809 12.461091995239258
11.428731918334961 10.658792495727539 9.816112518310547 8.862934112548828
7.930225849151611 6.979014873504639 5.906191825866699 4.732023239135742
3.70619797706604 2.8906478881835938 2.273542642593384 1.8709700107574463
1.5515344142913818 1.2728630304336548 1.0218982696533203
0.7814406156539917 0.595293402671814 0.4671807587146759 0.3404085636138916
0.29420292377471924 0.3925720453262329 0.4885733127593994
0.4826304614543915 0.4230322539806366 0.2793792486190796
0.0003629326820373535 -0.2482202649116516 -0.36572265625 -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- -- -- --],
mask = [False False False False False False False False False False False False
False False False False False False False False False False False False
False False False False False False False False False False False False
False True True False False True True True True True True True
True False False False False False True True True True True True
True True True True True True True True True True True True
True True True True True True True True True True True True
False False False False False False False False False False False False
False False False False False False False False False False False False
False False False False False False False False False False False False
False False False False False False False False False False False False
False False False False False False False False False False False False
False False False False False False False False False False False False
False False False False True True True True True True True True
True True True True True True True True True True True True],
fill_value = -1e+30)
確かに北極に近いところに-1000が入っている。sstはデータが格納されているdataだけでなく，True, Falseが格納されているmaskから構成されている。このような配列をMasked Arrayという。-1000もマスクすれば絵が描けそうだ。
In [5]:
import numpy.ma as ma
import netCDF4 as nc4
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
dataset = nc4.Dataset('HadISST1_SST_update.nc')
ice_covered = -1000
sst = ma.masked_values(dataset.variables['sst'][0, :, :], ice_covered)
lat = dataset.variables['latitude'][:]
lon = dataset.variables['longitude'][:]
ax = plt.axes(projection=ccrs.PlateCarree())
plt.contourf(lon, lat, sst, 60, transform=ccrs.PlateCarree(), vmin = -10.0, vmax = 35.0)
plt.colorbar()
ax.coastlines()
plt.show()
テキストデータ¶
NetCDFは軸の情報などが付加されていて便利であるが，HadISSTにはテキストデータで提供されているものがある。2004年のデータを読んでみよう。取得し解凍してから，1行読んでみる。
In [6]:
fin = open('HadISST1_SST_2004.txt', 'rt')
fin.readline()
Out[6]:
' 1 1 2004 180 rows 360 columns\n'
データ形式のとおり，最初の行はヘッダのようだ。次の行は北極に一番近い緯度円（89.5N）のデータ。
In [7]:
fin.readline()
Out[7]:
' -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000\n'
全て海氷だ。少し読み飛ばしてみる。
In [8]:
for j in range(20):
fin.readline()
fin.readline()
Out[8]:
'-32768-32768 -180 -180 -180 -180 -180 -180 -180 -180 -180 -180 -132 -58-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768 -180 -180 -180 -180 -180 -180 -180 -180-32768-32768 -180 -180 -180 -180 -180 -180 -180-32768-32768 -180 -180-32768-32768-32768-32768-32768-32768 -180 -1000-32768-32768-32768-32768 -180 -180 -180 -1000 -1000 -1000 -180 -1000-32768-32768-32768-32768-32768-32768 -180 -180 -1000 -1000 -1000 -180 -180 -180 -180 -180 -167 -64 15 34 38 18-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768 -1000 -1000 -180 59 82 149 211 247 259 261 268 275 275 264 245 217 186 185 213 248 289 332 373 417 464 515 569 592 587 592 604 623 645 657 662 665 668 670 666 654 633 609 591 578-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768 235 210 189 174 161 113 53 63 77 76 74 9 68 3-32768-32768 -141 -125-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768 -151 -134-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768-32768\n'
-32768は陸地で海面水温のないところ。NetCDFではmaskがTrueとなっていた。面倒なことに-32768の前には空白が省略されている。-1000に置き換えることにする。テキストの置き換えには，標準モジュールreを使う。
In [9]:
import re
land = -32768
seaice = -1000
re.sub(str(land), ' '+str(seaice), fin.readline())
Out[9]:
' -1000 -1000 -1000 -1000 -1000 -180 -180 -180 -180 -180 -180 -180 -170 -98 -58 -11 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -180 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -180 -180 -1000 -1000 -1000 -1000 -180 -180 -1000 -180 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -180 -180 -180 -160 -123 -12 43 59 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -180 -180 -180 -117 -59 -16 93 205 220 303 350 366 369 366 364 356 347 339 320 293 254 210 202 232 271 317 366 412 455 494 540 592 619 623 632 647 665 687 701 705 705 702 696 679 654 627 598 583 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 104 94 67 -1000 -7 -7 -27 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -133 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000 -1000\n'
緯度円上のデータが一つの文字列になっている。空白で区切り文字列をリストに分解する。
In [10]:
re.sub(str(land), ' '+str(seaice), fin.readline()).split()
Out[10]:
['-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-180', '-180', '-161', '-170', '-115', '-50', '-82', '-180', '-180', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-180', '-180', '-1000', '-1000', '-1000', '-180', '-180', '-180', '-180', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-180', '-180', '-137', '-81', '-14', '44', '82', '90', '99', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-180', '97', '1', '150', '254', '226', '248', '306', '397', '451', '459', '465', '421', '430', '428', '414', '397', '375', '356', '341', '323', '344', '319', '285', '270', '278', '301', '335', '382', '434', '480', '517', '556', '599', '626', '640', '659', '682', '705', '728', '743', '748', '742', '725', '707', '684', '656', '627', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '41', '-1000', '-1000', '-1000', '-1000', '-1000', '45', '65', '43', '21', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000', '-1000']
リストをnumpyの整数の配列に型変換する。
In [11]:
import numpy as np
np.array(re.sub(str(land), ' '+str(seaice), fin.readline()).split()).astype(int)
Out[11]:
array([ -62, 13, 88, -1000, -1000, -1000, -1000, -1000, -82,
-90, -115, -90, -132, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -166, -180, -180,
-1000, -1000, -180, -180, -180, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -180, -180, -1000, -1000,
-180, -180, -180, -108, -17, 57, 110, 135, 134,
124, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, 193, 293, 343,
357, 282, 251, 315, 387, 437, 459, 488, 537,
581, 589, 563, 536, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, 447, 430, 404, 385, 381,
386, 395, 423, 467, 506, 538, 571, 607, 638,
665, 692, 716, 737, 757, 768, 769, 757, 734,
711, 687, 659, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, 209, 197, 64, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -8,
-4, 90, 87, 85, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000,
-1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000, -1000])
以上をまとめて関数にする。
In [12]:
import re
import numpy as np
import numpy.ma as ma
nx = 360; ny = 180
def read_hadisst(fname, m):
fin = open(fname, 'rt')
sst = ma.zeros((m, ny, nx))
land = -32768
seaice = -1000
for k in range(m):
fin.readline()
for j in range(ny):
line = re.sub(str(land), ' '+str(seaice), fin.readline())
sst[k, j, :] = ma.masked_equal(np.array(line.split()).astype(int), seaice) * 0.01
fin.close()
return sst
In [13]:
sst = read_hadisst('HadISST1_SST_2004.txt', 12)
In [14]:
sst.shape
Out[14]:
(12, 180, 360)
In [15]:
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
ax = plt.axes(projection=ccrs.PlateCarree())
plt.contourf(lon, lat, sst[0,:,:], 60, transform=ccrs.PlateCarree(), vmin = -10.0, vmax = 35.0)
plt.colorbar()
ax.coastlines()
plt.show()
