気象庁MGDSST

気象庁の海面水温データ（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