用NCL处理卫星数据的教程【转】

mofangbao · 发表于 2012-2-9 08:58:01

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

转载自 figureacg
最终编辑 zhjstef
=====================
这个教程主要介绍如何用NCL读取卫星数据，如何将数据存储为二进制格式，如何作图。

这个数据是AMSR-E的L2A产品，更多信息请查阅相关文档。

在开始这个教程前，请确认已经安装好了NCL。

1. 运行NCL
和GrADS类似，NCL也可以运行在交互式环境中，即输入一条命令执行一条命令。
交互式环境使用：
在命令行下输入
$ ncl

退出请输入命令
ncl 0> exit

交互式环境不方便重复操作，一旦退出程序就需要全部重新来过，所以推荐编写NCL脚本。
一个NCL脚本主要包括以下几个部分：

load  外部函数

begin

;命令和操作

end

load用来载入外部函数。
文件操作、数据操作、数据可视化等代码都放在load 语句之后。
半角分号“;”后的内容是注释，NCL只有行注释功能。
begin 和 end 可有可无，但如果出现了begin，就必须有end 与之对应。

NCL脚本扩展名为.ncl。比如有一NCL脚本文件为 script.ncl，运行时可以用下面的方法：
$ ncl < script.ncl
或 $ ncl script.ncl

下面来写一个简单的NCL脚本。打开vi之类的编辑软件，将下面的内容复制进去，然后保存到卫星数据文件所在目录，文件名为test.ncl

;********************************************************
;  Load NCL scripts
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"
; load语句载入了NCL常用的几个外部程序，建议以后自己写脚本时也保留这6行的内容。
;********************************************************

  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  print(fi)

  end

运行这个脚本
$ ncl test.ncl > log

2. 熟悉数据
打开log文件
从第11行相关信息可以知道该文件为 HDFEOS_V2.5
第663行开始可以得知该卫星数据的生成时间等信息
第957行开始是卫星数据的维度信息(dimensions)
第983行起是变量信息(variables)

这些对文件内容进行描述的信息被称作元数据(metadata)，元数据对于数据的存储和发布很重要。它记录了数据集的主要特征，包括数据结构和内容描述，以及数据集开发的历史等。

3. 读取数据
下面以89V的数据为例
89V的元数据为：
short 89_0V_Res_1_TB__4 ( DataTrack_lo_Low_Res_Swath, DataXTrack_lo_Low_Res_Swath )
SCALE_FACTOR : 0.01
OFFSET : 327.68
UNIT : Kelvin
hdf_name : 89.0V_Res.1_TB
hdf_group : Low_Res_Swath/Data Fields
hdf_group_id : 4

short：该数据是short格式
89_0V_Res_1_TB__4：变量名
DataTrack_lo_Low_Res_Swath
DataXTrack_lo_Low_Res_Swath：数据大小，可以从维度信息(dimensions)了解数据大小。
SCALE_FACTOR
OFFSET：为了减小存储数据所占的空间，一般都会对实型数据进行加减乘除处理，以整型来存储(pack)，用户在使用时，需要将整型转换为实型(unpack)，会用到scale_factor和offset信息。

修改begin 和 end 之间的内容：
  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  data =new(dimsizes(fi->89_0V_Res_1_TB__4),float,-9999)
  data = (fi->89_0V_Res_1_TB__4)* 0.01 + 327.68
  printMinMax(data,False)
  end

这样，data变量里存储的就是实型的89V数据。printMinMax输出data的最大值和最小值。
运行这个脚本，显示结果为

4. 作图
现在我们有89V数据了，但是还没有经纬度信息。
从元数据中发现多组经纬度信息：
Latitude__3，Latitude__7，Latitude__11
Longitude__3，Longitude__7，Longitude__11
从各个变量的数据大小得知应该选Latitude__3和Longitude__3
89_0V_Res_1_TB__4 ( DataTrack_lo_Low_Res_Swath, DataXTrack_lo_Low_Res_Swath )
Latitude__3 ( DataTrack_lo_Low_Res_Swath, DataXTrack_lo_Low_Res_Swath )
Longitude__3 ( DataTrack_lo_Low_Res_Swath, DataXTrack_lo_Low_Res_Swath )

将脚本中begin和end中间的部分改为：

  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  data =new(dimsizes(fi->89_0V_Res_1_TB__4),float,-9999)
  data = (fi->89_0V_Res_1_TB__4)* 0.01 + 327.68

  lat2d = fi->Latitude__3
  lat2d@units ="degrees_north"
  lon2d = fi->Longitude__3
  lon2d@units ="degrees_east"

  data@lat2d = lat2d
  data@lon2d = lon2d

  end

这样就将数据和经纬度信息对应起来。

将begin和end之间的内容改为
  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  data =new(dimsizes(fi->89_0V_Res_1_TB__4),float,-9999)
  data = (fi->89_0V_Res_1_TB__4)* 0.01 + 327.68

  lat2d = fi->Latitude__3
  lat2d@units ="degrees_north"
  lon2d = fi->Longitude__3
  lon2d@units ="degrees_east"

  data@lat2d = lat2d
  data@lon2d = lon2d

wks = gsn_open_wks("ps" ,"AMSR_E")
gsn_define_colormap(wks,"WhBlGrYeRe")

  res                      = True

  res@mpFillOn             = False
  res@cnFillOn             = True
  res@cnFillMode          = "RasterFill"
  res@cnLinesOn          = False
  res@cnLineLabelsOn       = False
  res@cnLevelSelectionMode  = "ManualLevels"
  res@cnMinLevelValF       = 160
  res@cnMaxLevelValF       = 300
  res@cnLevelSpacingF    = 5
  res@gsnSpreadColors    = True
  res@lbLabelBarOn       = True
  res@lbLabelAutoStride    = True
  res@lbBoxLinesOn       = False

  plot  = gsn_csm_contour_map(wks,data,res)

  end

运行脚本(这个过程比较慢)，当前目录下会生成AMSR_E.ps文件。
用GS软件打开

看来这个地图投影方式不太合适，换一个投影方式试试。
将begin和end之间的内容改为：

  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  data =new(dimsizes(fi->89_0V_Res_1_TB__4),float,-9999)
  data = (fi->89_0V_Res_1_TB__4)* 0.01 + 327.68

  lat2d = fi->Latitude__3
  lat2d@units ="degrees_north"
  lon2d = fi->Longitude__3
  lon2d@units ="degrees_east"

  data@lat2d = lat2d
  data@lon2d = lon2d

wks = gsn_open_wks("ps" ,"AMSR_E")
gsn_define_colormap(wks,"WhBlGrYeRe")

  res                      = True

  res@mpFillOn             = False
  res@cnFillOn             = True
  res@cnFillMode          = "RasterFill"
  res@cnLinesOn          = False
  res@cnLineLabelsOn       = False
  res@cnLevelSelectionMode  = "ManualLevels"
  res@cnMinLevelValF       = 160
  res@cnMaxLevelValF       = 300
  res@cnLevelSpacingF    = 5
  res@gsnSpreadColors    = True
  res@lbLabelBarOn       = True
  res@lbLabelAutoStride    = True
  res@lbBoxLinesOn       = False

  res@mpProjection             = "Satellite" ; choose map projection
  res@mpCenterLonF             = 130.0    ; choose center lon
  res@mpCenterLatF             = 25.       ; choose center lat
  res@mpSatelliteDistF          = 5.0       ; choose satellite view

  plot  = gsn_csm_contour_map(wks,data,res)

  end

运行脚本，然后再看看。

5. 存储数据
存储为FORTRAN unformatted direct 文件：

  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  data =new(dimsizes(fi->89_0V_Res_1_TB__4),float,-9999)
  data = (fi->89_0V_Res_1_TB__4)* 0.01 + 327.68

  lat2d = fi->Latitude__3
  lon2d = fi->Longitude__3

  fo1 = "./89V.dat"
  fbindirwrite(fo1,data )
  fo2 = "./lat.dat"
  fbindirwrite(fo2,lat2d )
  fo3 = "./lon.dat"
  fbindirwrite(fo3,lon2d )

  end

当前目录生成89V.dat，lon.dat和lat.dat三个文件。读数据用的FORTRAN程序：
open(..., form="unformatted", access="direct", recl= )
------------------------------
存储为netCDF格式：

  begin

  fi = addfile("./AMSR_E_L2A.hdf","r")
  data =new(dimsizes(fi->89_0V_Res_1_TB__4),float,-9999)
  data = (fi->89_0V_Res_1_TB__4)* 0.01 + 327.68

  lat2d = fi->Latitude__3
  lon2d = fi->Longitude__3

  fo = addfile("./data.nc","c")
  fo->89V = data
  fo->lat2d = lat2d
  fo->lon2d = lon2d

  end

当前目录生成data.nc文件
可以用ncdump等工具看文件内容。

mofangbao · 发表于 2012-2-19 23:23:55

灰太狼发表于 2012-2-19 11:28

登录/注册后可看大图

mofangbao：您好！请问按照您的方法用NCL可以处理BUFR格式的卫星亮温数据吗？如gdas1.t00z.1bamua.tm00.buf ...

抱歉啦，这个是转载的帖子，我并没有使用过，你可以自己试一下的

freekiller · 发表于 2012-2-11 11:09:41

学习了，多谢

zhaoyanlainj · 发表于 2012-2-18 10:13:37

多谢，学习啦

qxtlyf · 发表于 2012-2-19 10:30:59

好东西，谢谢了

灰太狼 · 发表于 2012-2-19 11:28:02

mofangbao：您好！请问按照您的方法用NCL可以处理BUFR格式的卫星亮温数据吗？如gdas1.t00z.1bamua.tm00.bufr_d

zhf1268 · 发表于 2012-2-20 15:42:18

呵呵，学习学习！

liushuyanhe · 发表于 2012-7-6 02:22:54

非常感谢。学习啦

qxtlyf · 发表于 2012-7-6 07:18:14

呵呵，学习学习

bapibenben · 发表于 2012-7-7 12:54:18

帮助学习ncl

		自动登录	找回密码
密码			立即注册

[经验总结] 用NCL处理卫星数据的教程【转】

登录后查看更多精彩内容~