気象庁の海面水温データ（0.25x0.25度、daily)であるMGDSSTをウェブからダンロードしてから、データを読み込むまでのモジュールをつくる。MGDSSTはNEAR-GOOSのサイトから入手が可能である。
idとpasswordが必要なので、あらかじめ登録しておく。
以下がスクリプト
usrname="xxxxx"
password="*****"
のところを自分のidとpasswordで置き換える。
import numpy as np
import os
def downloadfile(infile="remotefile",outfile="localfile"):
"""
download file from NEAR-GOOS data site with id and password
input
infile: remote file on the web
outfile: local file
set password and id below
"""
import urllib2
host="goos.kishou.go.jp:80"
realm="NEAR-GOOS"
usrname="xxxxx"
password="*****"
auth_handler = urllib2.HTTPBasicAuthHandler()
auth_handler.add_password(realm, host, usrname, password)
opener = urllib2.build_opener(auth_handler)
urllib2.install_opener(opener)
f=urllib2.urlopen(infile)
open(outfile, 'w').write(f.read())
#----------------------------------------------------------
def readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)):
"""
load MGDSST data from NEAR-GOOSE site
read data with fortran (f2py)
select region specified by lonrange and latrange
Usage:
lon,lat,sst= readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)
Input
dir: directory for data download
year,month,day: date of the data
lonrange,latrange = longitude and latitude range for data
output
lon, lat: longitude, latitude
sst: sst data
"""
# date
import datetime
t=datetime.datetime(year,month,day)
syear=t.strftime("%Y")
smon=t.strftime("%b")
sday=t.strftime("%d")
# filename
dataname="mgdsst."+smon+sday
filename=dir+"/"+dataname+"."+syear
filenamegz=filename+".gz"
print filename
urlname="http://goos.kishou.go.jp/rrtdb/usr/pub/JMA/mgdsst/"+syear+"/"+dataname+".gz"
# file exist?
if not os.path.exists(filename):
if not os.path.exists(filenamegz):
downloadfile(infile=urlname,outfile=filenamegz)
os.system("gunzip %(filenamegz)s" %locals())
# temp.bin is used for temporary file
if os.path.lexists("temp.bin"):
raise IOError("temporary file temp.bin somehow exist! Remove it!")
else:
os.system("ln -s %(filename)s temp.bin" %locals())
# data read through fortran (f2py)
from read_sst import read_sst
im=1440
jm=720
sst=np.zeros([jm,im])
sst,ierr=read_sst()
assert ierr==0, "data error in fortran"
lon=0.125+0.25*np.arange(im)
lat=-89.875+0.25*np.arange(jm)
os.system("rm temp.bin")
# range selection
xrange=np.logical_and(lon>=lonrange[0],lon<=lonrange[1])
yrange=np.logical_and(lat>=latrange[0],lat<=latrange[1])
sst=(sst[yrange,:])[:,xrange]
return lon[xrange],lat[yrange], sst
#-----------------------------------------------------------------------
if __name__ == '__main__':
import cdms2,MV2,vcs #cdat
# read data
lon,lat,sst=readMGDSST(dir="DATA",year=2011,month=1,day=1,lonrange=(0,360),latrange=(-90,90))
# hearafter plot with cdat
# First let's create a variable
var = MV2.masked_array(sst,mask=(sst<-990.),id="SST")
# Now create the axis
lons=cdms2.createAxis(lon)
lats=cdms2.createAxis(lat)
# And decorate it
lons.id = 'longitude'
lats.id = 'latitude'
lons.units = 'degrees_east'
lats.units = 'degrees_north'
# Just to make sure we will describe this axis as a longitude one
lons.designateLongitude()
lats.designateLatitude()
# Now let's apply these to the variable
var.setAxisList((lats,lons))
#var.info() to see info
#plot
v = vcs.init()
v.plot(var)
v.png("MGDSST.png")
del v
上のサンプルメインプログラムでは、2011年1月1日のデータをディレクトリ"DATA"にダウンロードし、経度0から360度、緯度-90度から90度の範囲で読み込んでいる。
その上で、cdatを用いて可視化した。
結果は以下の図。
ダウンロードしたバイナリファイルを読むのにfortranのプログラムをモジュールとしてつかう。清水慎吾さんのデータを読むサンプルプログラムを参考に、以下のようなfortranサブルーチンを作った（配列はpython風に緯度と経度を入れ替えている。欠損値は-999.）。
このプログラムをf2pyを使って以下のようにコンパイルすればpythonのモジュールとして使用が可能になる（上のプログラムのread_sst)。
f2py --f90exec=gfortran -c -m read_sst read_sst.f90
(ここではgfortranを使っている。）
read_sst.ｆ90
subroutine READ_SST(sst,ierr)
implicit none
integer :: year,mon,day
integer :: i,j,jj
integer,parameter :: im=1440
integer,parameter :: jm=720
!
real,intent(out):: sst(jm,im)
integer, intent(out):: ierr
!
integer :: isst(im,jm),ist
real :: undef = -999.
open(55,file="temp.bin",status='old',form='formatted',iostat=ierr,err=999)
read(55,*) year,mon,day
write(6,*) year,mon,day
do j=1,jm
read(55,10,iostat=ierr,err=999) (isst(i,j),i=1,im)
enddo
close(55)
10 format(1440i3)
do j=1,jm
do i=1,im
jj = jm -j +1
ist = isst(i,jj)
if(ist .eq. 888 .or. ist.eq.999) then
sst(j,i) = undef
else
sst(j,i) = real(ist)*0.1
endif
enddo
enddo
return
999 print *, "Read Error!"
end subroutine READ_SST
