Argoデータの利用

・フロートの情報、取得情報は全てファイルの格納されています.

・netCDFフォーマットになっています(拡張子は.nc).

Argoデータのダウンロード

本当はいろいろあるようですが、私もそんなに知っているわけではないです.

フロートの位置をみる. csvやjson形式でダウンロードすることも可能らしいです. 私はArgovisが使いやすかったです. (期間の指定、管理している研究機関がどこか見やすい): https://argo.ucsd.edu/data/data-visualizations/

データのダウンロード: https://data-argo.ifremer.fr/

ディレクトリの構造

• ダウンロード先にアクセスすると、色々なフォルダがあります

dac: 研究機関別にデータを保存

geo: 海域別 (atlantic, indian, pacific) にデータを保存

latest: 直近1ヶ月に新規作成・更新されたデータを保存

個々のフロートを調べたい場合、dacからデータを探すことが多い気がします. なので今回はとりあえず、dac/jma/2902474 という適当なファイルを開いてみます.

• ファイルの種類

2902474_Rtraj.nc: 軌跡のデータ

2902474_meta.nc: 投入情報・搭載センサの情報など

2902474_prof.nc: プロファイルデータのまとめ（こっちを使ってもいいような？）

2902474_tech.nc: バッテリーの残量など、技術的な情報

profiles/: 個々のプロファイルがいっぱい入ってる

import numpy as np
import xarray as xr  #netCDF4を使う人も多いです

nc = xr.open_dataset('D2902474_000.nc')
nc

<xarray.Dataset>
Dimensions:                       (N_PROF: 1, N_PARAM: 3, N_LEVELS: 116, N_CALIB: 2, N_HISTORY: 17)
Dimensions without coordinates: N_PROF, N_PARAM, N_LEVELS, N_CALIB, N_HISTORY
Data variables: (12/64)
    DATA_TYPE                     object b'Argo profile    '
    FORMAT_VERSION                object b'3.1 '
    HANDBOOK_VERSION              object b'1.2 '
    REFERENCE_DATE_TIME           object b'19500101000000'
    PLATFORM_NUMBER               (N_PROF) object b'2902474 '
    PROJECT_NAME                  (N_PROF) object b'J-ARGO                   ...
    ...                            ...
    HISTORY_ACTION                (N_HISTORY, N_PROF) object b'IP  ' ... b'    '
    HISTORY_PARAMETER             (N_HISTORY, N_PROF) object b'              ...
    HISTORY_START_PRES            (N_HISTORY, N_PROF) float32 nan nan ... nan
    HISTORY_STOP_PRES             (N_HISTORY, N_PROF) float32 nan nan ... nan
    HISTORY_PREVIOUS_VALUE        (N_HISTORY, N_PROF) float32 nan nan ... nan
    HISTORY_QCTEST                (N_HISTORY, N_PROF) object b'              ...
Attributes:
    title:                Argo float vertical profile
    institution:          JAMSTEC
    source:               Argo float
    history:              2014-05-08T11:20:30Z creation;2014-06-20T00:22:26Z ...
    references:           http://www.argodatamgt.org/Documentation
    comment:              free text
    user_manual_version:  3.1
    Conventions:          Argo-3.1 CF-1.6
    featureType:          trajectoryProfile

• Dimension

・N_PROF: プロファイルの数、だいたい1個

・N_PARAM: パラメータの数 (p,T,S)

・N_LEVELS: 深さ方向の層の数

• Data variables

多すぎるので大事なものだけ・・・

・PRES_ADJUSTED: 補正圧力

・TEMP_ADJUSTED: 補正水温

・PSAL_ADJUSTED: 補正塩分

・LATITUDE: 緯度

・LONGITUDE: 経度

・JULD: 時間(1950.01.01から)

QC (Quality Control): たまに変なデータが紛れ込んでいるので、その判定をしています.

・PRES_ADJUSTED_QC: 補正圧力

・TEMP_ADJUSTED_QC: 補正水温

・PSAL_ADJUSTED_QC: 補正塩分

・POSITION_QC: 緯度

・JULD_QC: 時間(1950.01.01から)

lat = nc['LATITUDE']
lon = nc['LONGITUDE']
date = nc['JULD']
p = nc['PRES_ADJUSTED'][:]
T = nc['TEMP_ADJUSTED'][:]
S = nc['PSAL_ADJUSTED'][:]

### 余計な次元は削減 ###
p = np.squeeze(p)
T = np.squeeze(T)
S = np.squeeze(S)

### 鉛直プロファイルを書いてみよう! ###
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(9,6))
ax = fig.add_subplot(1,2,1)
ax.plot(T,p,color='k')
ax.set_xlabel('In-situ temperature [$\circ$C]',fontsize=13)
ax.set_ylabel('Pressure [dBar]',fontsize=13)
ax.tick_params(axis='both',labelsize=13)
ax.grid(which='major',color='gray',linestyle='dashed')
ax.invert_yaxis()

ax = fig.add_subplot(1,2,2)
ax.plot(S,p,color='k')
ax.set_xlabel('Salinity [psu]',fontsize=13)
ax.set_ylabel('Pressure [dBar]',fontsize=13)
ax.tick_params(axis='both',labelsize=13)
ax.grid(which='major',color='gray',linestyle='dashed')
ax.invert_yaxis()

plt.tight_layout()
plt.show()

ポテンシャル密度を求める

海洋学で使いそうな物理量を計算するtoolbox: https://www.teos-10.org/software.htm

anacondaでインストールする

conda install -c conda-forge gsw

import gsw

### 絶対塩分 ###
SA = gsw.SA_from_SP(S,p,lon,lat)

### potential density ###
pd = gsw.pot_rho_t_exact(SA,T,p,0.0)-1000.0

### 鉛直プロファイル ###
fig = plt.figure(figsize=(9,6))
ax = fig.add_subplot(1,1,1)
ax.plot(pd,p,color='k')
ax.set_xlabel('Potential density [$kg/m^3$]',fontsize=13)
ax.set_ylabel('Pressure [dBar]',fontsize=13)
ax.tick_params(axis='both',labelsize=13)
ax.grid(which='major',color='gray',linestyle='dashed')
ax.invert_yaxis()
plt.show()

フロートを追跡する

1点だけ見ることはあまりないと思うので、フロートを追跡して断面図を描きたいですね.

今回は2902474_prof.ncを使います.

nc = xr.open_dataset('2902474_prof.nc')
nc

<xarray.Dataset>
Dimensions:                       (N_PROF: 146, N_PARAM: 3, N_LEVELS: 116, N_CALIB: 2, N_HISTORY: 0)
Dimensions without coordinates: N_PROF, N_PARAM, N_LEVELS, N_CALIB, N_HISTORY
Data variables: (12/64)
    DATA_TYPE                     object b'Argo profile    '
    FORMAT_VERSION                object b'3.1 '
    HANDBOOK_VERSION              object b'1.2 '
    REFERENCE_DATE_TIME           object b'19500101000000'
    DATE_CREATION                 object b'20140103223437'
    DATE_UPDATE                   object b'20190918001051'
    ...                            ...
    HISTORY_ACTION                (N_HISTORY, N_PROF) object 
    HISTORY_PARAMETER             (N_HISTORY, N_PROF) object 
    HISTORY_START_PRES            (N_HISTORY, N_PROF) float32 
    HISTORY_STOP_PRES             (N_HISTORY, N_PROF) float32 
    HISTORY_PREVIOUS_VALUE        (N_HISTORY, N_PROF) float32 
    HISTORY_QCTEST                (N_HISTORY, N_PROF) object 
Attributes:
    title:                Argo float vertical profile
    institution:          FR GDAC
    source:               Argo float
    history:              2019-09-18T00:10:51Z creation
    references:           http://www.argodatamgt.org/Documentation
    user_manual_version:  3.1
    Conventions:          Argo-3.1 CF-1.6
    featureType:          trajectoryProfile

lat = nc['LATITUDE']
lon = nc['LONGITUDE']
date = nc['JULD']
p = nc['PRES_ADJUSTED'][:]
T = nc['TEMP_ADJUSTED'][:]
S = nc['PSAL_ADJUSTED'][:]

### 余計な次元は削減 ###
p = np.squeeze(p)
T = np.squeeze(T)
S = np.squeeze(S)

流跡

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter

plt.figure(figsize=(10,8))
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.add_feature(cfeature.LAND)
ax.coastlines(lw=0.5)

ax.set_title('2902474',fontsize=24)

### 目盛を描く緯度経度の値を設定 ###
ax.set_xticks(np.arange(130,150.1,5),crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(24,40,4),crs=ccrs.PlateCarree())
    
### 目盛の表示形式を度数表記にする ### 
latfmt = LatitudeFormatter()
lonfmt = LongitudeFormatter(zero_direction_label=True)
ax.xaxis.set_major_formatter(lonfmt)
ax.yaxis.set_major_formatter(latfmt)
ax.axes.tick_params(labelsize=18)
ax.set_extent([130,150,24,38],ccrs.PlateCarree())
    
ax.plot(lon,lat,transform=ccrs.PlateCarree())
ax.plot(lon,lat,'+k',transform=ccrs.PlateCarree(),markersize=7)
ax.plot(lon[0],lat[0],'ok',transform=ccrs.PlateCarree(),markersize=7)
    
plt.show()    

鉛直断面を描く

各プロファイルで微妙に深度が違うので、補間によって深度を揃えます

### 不良データを取り除いて内挿する ###
from scipy import interpolate
T_itp = []
prs_ref = np.arange(10,2010,10)

for j in range(0,146):
    temp = T[(j)]
    pres = p[(j)]
    dat0 = np.vstack([pres,temp])
    dat = dat0[:,~np.isnan(dat0).any(axis=0)]
    
    f = interpolate.Akima1DInterpolator(dat[0,:],dat[1,:])
    T_itp.append(f(prs_ref))

T_itp = np.squeeze(np.array(T_itp))

混合層深度

ポテンシャル密度や水温で定義します. 今回は水温が10mから0.5℃低下した深度とします.

mld = np.zeros([T.shape[0]])
for i in range(0,T.shape[0]):
    norm = np.array(np.where(T_itp[i,:] < T_itp[i,0]-0.5))[0,0]
    mld[i] = prs_ref[norm]

### gridの作成 ###
Y,X = np.meshgrid(prs_ref,date)

fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(1,1,1)
ax.set_title('2902474',fontsize=24)

ax.plot(date,mld,color='black',label='MLD',linestyle='dashed')
sc = ax.pcolormesh(X,Y,T_itp,cmap='bwr',vmin=0,vmax=30,shading='nearest')
ax.set_ylabel('Pressure [dBar]',fontsize=18)
ax.tick_params(axis='both',labelsize=18)
ax.invert_yaxis()
       
cb = plt.colorbar(sc, pad=0.15, aspect=35, shrink=0.8,orientation="horizontal")
cb.set_label("[$^\circ$C]",labelpad=5,size=18)
cb.ax.tick_params(labelsize=18)

plt.show()