求助：处理 fwrite提取 wrf输出的dat文件 中的单要素二进制文件的问题

zhudan0723

我的具体情况是这样：
首先，从wrf已经生成的wrfout_d02.dat，用Grads的fwrite 提取 tk的一个层次如（850），一个时次的二进制数据文件。
ctl 文件内容：
dset K:\WRF-OUT\2012050912\3km\wrfout_d02.dat
undef 1.e35
pdef 312 285 lcc  39.511  121.599    1.000    1.000  60.  30.  126.500   3000.   3000.
xdef 1098 linear  119.9   0.01351351
ydef  748 linear   38.4   0.01351351
zdef    7 levels  
1000.00000
  925.00000
  850.00000
  700.00000
  500.00000
  300.00000
  200.00000
tdef        85 linear 12z09may2012  1hr
vars   27
UMET             7 0 U Compoment of wind - rotated (diagnostic)        
VMET             7 0 V Component of wind - rotated (diagnostic)        
W                7 0 W Component of wind                              
THETA            7 0 Theta                                             
TK               7 0 Temperature in K                                 
Z                7 0 Height (m)                                       
QVAPOR           7 0 Vapor (kg/kg)                                    
TD               7 0 Dewpoint Temperature in C (diagnostic)            
RH               7 0 Relative Humidity (diagnostic)                    
RAINC           0  0 ACCUMULATED TOTAL CUMULUS PRECIPITATION           
RAINNC          0  0 ACCUMULATED TOTAL GRID SCALE PRECIPITATION        
slvl            0  0 sea level pressure                                
T2              0  0 TEMP at 2 M                                       
U10M            0  0 U at 10 M - rotated                              
V10M            0  0 V at 10 M - rotated                              
XLAT            0  0 LATITUDE, SOUTH IS NEGATIVE                       
XLONG           0  0 LONGITUDE, WEST IS NEGATIVE                       
PSFC            0  0 SURFACE PRESSURE(PA)                              
RH2             0  0 Relative humidity(%) at 2m (diagnostic)           
TD2             0  0 Dewpoint temperature(K) at 2m (diagnostic)        
T2MAX           0  0 Maximum temperature(K) at 2m (diagnostic)         
T2MIN           0  0 Minimum temperature(K) at 2m (diagnostic)         
HIGHCLD         0  0 High cloud fraction[%] (diagnostic)               
MIDCLD          0  0 Middle cloud fraction[%] (diagnostic)            
LOWCLD          0  0 Low cloud fraction[%] (diagnostic)               
TOTCLD          0  0 Total cloud fraction[%] (diagnostic)              
VIS             0  0 Visibility in km (diagnostic)                     
endvars
提取二进制文件的的gs
'c'
'set lev 850'
'set t 25'
'set gxout fwrite'
'set fwrite D:\wrfGRD\GrdOut\TK\850\2012050920.025.grd'
'd TK'
'disable fwrite'
===============================
然后，我想通过程序，通过一个站点的经纬度，找到一个与它最接近的个点，认作是该站点的数据（精度为3km的，我认为可以近似，就不想插值了）
根据，
xdef 1098 linear  119.9   0.01351351
ydef  748 linear   38.4   0.01351351
在程序中，设计的个点坐标为：X = 119.9 + (i - 1) * 0.01351351   Y = 38.4 + (j - 1) * 0.01351351 （i[1,1098]，j[1,748]是两个循环变量）
最后，我读出的点 温度为10.6956摄氏度，站点坐标 125.22432294，43.89999857
格点坐标  134.7243204          48.49459197
但是在原始图上,看大约是8度。
在网上查，说是wrfout_d02.ctl有投影，是投影的原因。
求论坛中的高手详解！

zhudan0723

各位高手，请帮帮忙，教教在下~

catmiaow

zhudan0723 发表于 2012-5-11 13:35

                               
登录/注册后可看大图

各位高手，请帮帮忙，教教在下~

直接在GS里加站点的经纬度就可以了吧。
'set lev 850'
'set t 25'
'set lon ****'
'set lat ***'
'set gxout fwrite'
'set fwrite D:\wrfGRD\GrdOut\TK\850\2012050920.025.grd'
'd TK'
'disable fwrite'

zhudan0723

catmiaow 发表于 2012-5-11 20:35

                               
登录/注册后可看大图

直接在GS里加站点的经纬度就可以了吧。
'set lev 850'
'set t 25'

我试试！不过设置经纬度，那生成的二进制文件，经纬度起始坐标和间隔怎么知道！？

我看到动力论坛里有说
set x xxxx
set y xxxx

天目神眉

最后，我读出的点 温度为10.6956摄氏度，站点坐标 125.22432294，43.89999857
格点坐标  134.7243204          48.49459197
但是在原始图上,看大约是8度。
没太明白楼主的意思，你算的这个站点坐标和格点坐标相差就很大，应该不止三公里吧？另外，原图上也不是8度啊。
你可以查看一下我前期的帖子，有这个的具体过程，只是和你的方法不太一样！可以参考一下！

zhudan0723

天目神眉 发表于 2012-5-11 22:18

                               
登录/注册后可看大图

最后，我读出的点 温度为10.6956摄氏度，站点坐标 125.22432294，43.89999857
格点坐标  134.7243204     ...

不好意思！贴的时候没注意，错了！
格点坐标是根据ctl算的终止格点坐标
而站点坐标是格点坐标。站点坐标为 125.22 43.9

catmiaow

zhudan0723 发表于 2012-5-11 22:10

                               
登录/注册后可看大图

我试试！不过设置经纬度，那生成的二进制文件，经纬度起始坐标和间隔怎么知道！？

我看到动力论坛里有 ...

经纬度起始坐标和间隔ctl里就有，我以前提取过，用着很笨的方法，先是用grads把数据转化一下，主要是因为WRF生成的数据有pdef,而且存放顺序为反序，直接用Fortran提取就会出现读入错误（也可能是我的Fortran没编好）。以下是我用的GS：
'reinit'
'open G:\exercise\ss\tc.ctl'
'set gxout fwrite'
'set fwrite G:\exercise\ss\test_16.dat'
tt=1
while(tt<=169)
'set t ' tt
hh=1
while(hh<=19)
'set z 'hh
'd tc'
hh=hh+1
endwhile
tt=tt+1
endwhile
'disable fwrite'
'reinit'

Fortran提取数据：
program Main
implicit none
integer, parameter :: iLon = 179, iLat = 132, iDay = 169,zlev=19
real rHgt(iLon, iLat, zlev,iDay)
integer iRec
integer iX, iY, iT,iz

open (10, file = "test_16.dat", form = "unformatted", access = "direct", recl = iLon * iLat)
iRec = 0
do iT  = 1, iDay
   do iz=1,zlev
   iRec = iRec + 1
   read (10, rec = iRec) ((rHgt(iX, iY, iz,iT), iX = 1, iLon), iY = 1, iLat)
end do
end do
close (10)
do iT  = 1, iDay
do iz=1,zlev
do iY  =1,iLat
do iX  =1,iLon
if (rHgt(iX, iY,iz, iT)>100000)then
rHgt(iX, iY, iz,iT)=0
end if
end do
end do
end do
end do
open (20, file = "G:\exercise\ss\tc\tc16.txt")
do iT = 1, iDay
   do iz=1,zlev
      write (20,'(i3,f10.2)') iz,rHgt(63, 56, iz,iT)
end do
end do
close (20)

rHgt(63, 56, iz,iT)是所选择的网格点上的数据，63,56是在grads里用set x/y和set lat/lon来尝试寻找最接近我选的那个站点的网格点。
我这方法是很笨了。

zhudan0723

catmiaow ,感谢您~
“63,56是在grads里用set x/y和set lat/lon来尝试寻找最接近我选的那个站点的网格点”您的意思63 56这点是您自己看出来的~而 “iLon = 179, iLat = 132”是个点说，也就是说个点数和个点坐在的经纬度很难换算。是这样么？
“主要是因为WRF生成的数据有pdef,而且存放顺序为反序”这里能详细解释一下么~！？

zhudan0723

fwrite输出的数据个点数据
应该是按着图画出来的，
xdef 1098 linear  119.9   0.01351351
ydef  748 linear   38.4   0.01351351
算了一下单个数据文件里有823152个数据，我用的是单精度浮点型，vb中是Single，2字节
我试验一下 （823152-1）/1099=749
说明fwrite输出的数据个点数据网格是1099X749的，较ctl多1个，或者我理解错了。但是这样还是较总个点数多一个，不知道。
在ctl中用两个变量，xlat 和xlong的输出，我查了是经纬度，通过fwrite输出的这变量，找到的占地数据正确。
为什么网格会窜位，怎么多1个，和总网格数多出一个来，弄不明白了！
不过，现在利用xlat 和xlong这两变量可以找到我想要的了！如果转成micaps，格式也不成问题了！
看家园中有国家气象中心邓莲堂开发的一个WRF模式的后处理程序-可生成GrADS和Micaps格式 ，我对fortran了解不多，简单的程序还能看明白，这位高手的代码初略看了一下，感觉很复杂，语法内容也多，没看懂。也许，对症的~可惜我看不懂~  

catmiaow

zhudan0723 发表于 2012-5-12 10:36

                               
登录/注册后可看大图

catmiaow ,感谢您~
“63,56是在grads里用set x/y和set lat/lon来尝试寻找最接近我选的那个站点的网格点”您 ...

格点数与格点坐标不难算啊，我只是懒得算，你上面的公式是对的，但不知道为什么你算出来的点不对。63,56是根据后处理给出的CTL里面的起始经纬度和网格距来算，不是根据namelist里的网格点数算的。
ARWpost生成的ctl里有options  byteswapped，数据存放顺序取反序，还有这个一句pdef   69  57 lcc  30.579  103.923   35.000   29.000  60.00000  30.70000  103.97000   3000.000   3000.000，69+1,57+1是运行模式时namelist里的网格点数，30.579  103.923是中心经纬度，这语句具体什么意思我也没明白，总之要用ARWpost生成的数据，就根据ctl里的xdef  179 linear  102.71060   0.01351351
ydef  132 linear   29.68310   0.01351351来算。以前有为同学是用了namelist里的网格点算，结果不对，不知道你是不是这种情况。我知道的也就这些，能力有限。