画锋生函数的问题：直接用nc文件运算和改写成grd文件画出的结果不一样？

echo要读书

================06.09 补充一下问题========================
补充一下我想达到的效果以及方法。
我是想画出15年的特定时段（比如6-7月）平均锋生强度（根据文献，我将格点平均锋生强度叩定义为该格点上锋生函数为正的数值的累计值与锋生函数为正值事件的总频次的比值。），于是我在将nc转写grd文件过程中，先求出某年某时间段的平均锋生强度（但是在用if语句筛选的过程中貌似也出现了问题：if语句没有用？），再将15年的锋生函数平均。

由于我所需的时间年份跨度大（1991-2005），不能直接在gs文件里直接打开60多个文件，貌似我就必须得把所需的时间段的数据提到grd运算？同时又由于每年我所需的时间段不规则，所以我想直接在转换文件里面把平均锋生强度提出来写进grd文件。

故提取平均锋生强度的程序如下，主要出现的问题有两个：1.直接用nc进行运算的时候就没问题，改用grd的时候就会出现LZ那种问题；2.我用nc运算的时候，貌似if语句的筛选没有作用？画出来的图还是有负值的？


'set gxout fwrite'
'set fwrite Z:\prog\grd\fro\try.grd'

**********************************1991
'sdfopen z:\Guo\climate\temdaily\air.1991.nc'
'sdfopen z:\Guo\climate\rhumdaily\rhum.1991.nc'
'sdfopen Z:\Guo\climate\uwnd\uwnd.1991.nc'
'sdfopen Z:\Guo\climate\vwnd\vwnd.1991.nc'
'define all=0.0'

'set x 33 55'
'set y 41 59'
'set z 3'

tt=159
'define m=0.0'
'define day=0'
while (tt<=165)
'set t 'tt''
'define tem=air.1'
'define rh=rhum.2'
'define wu=uwnd.3'
'define wv=vwnd.4'

*****pot**********
'define prs=850'
'define es=(6.112*exp(17.67*(tem-273.15)/(tem-38)))'
'define q=rh*(0.62197*es/(prs-0.378*es))/100'
'define pot=tem*pow((1000/prs),(0.2854*(1.0-0.28*q)))'

****gradient***

'define r=6.37e6'
'define dtx = cdiff(pot,x)'
'define dty = cdiff(pot,y)'
'define dx = cdiff(lon,x)*r*cos(lat*3.1416/180)*3.1416/180'
'define dy = cdiff(lat,y)*r*3.1416/180'
'define ga=1/(mag(dtx/dx,dty/dy))'

******u,v winds************
'define dux=cdiff(wu,x)'
'define duy=cdiff(wu,y)'
'define dvx=cdiff(wv,x)'
'define dvy=cdiff(wv,y)'

'fr=-1*ga*((dtx/dx)*((dux/dx)*(dtx/dx)+(duy/dy)*(dty/dy))+(dty/dy)*((dvx/dx)*(dtx/dx)+(dvy/dy)*(dty/dy)))'

if(fr>=0)      ****此处if语句貌似没起作用？？？
'm=m+fr'
'day=day+1'
endif
tt=tt+1
endwhile

'd m/day'

;
==============================

我想求取一段时间的锋生函数的平均值，直接用nc文件画图得到量极为10e-10，这个结果应该是正确的。但是我将所需时间段的函数计算后单独写在一个grd文件中，画出的图就不对了（画出的图的数据都是-9e10+8），不知道我是哪里写错了呢？。。

直接用nc读写的gs文件：
'sdfopen z:\Guo\climate\pottem\pottmp.sig995.1991.nc'
'sdfopen Z:\Guo\climate\surfaceu\uwnd.sig995.1991.nc'
'sdfopen Z:\Guo\climate\surfacev\vwnd.sig995.1991.nc'

'set x 33 55'
'set y 41 59'

tt=225
'define m=0.0'
while (tt<=239)
'set t 'tt''
'define pot=pottmp.1'
'define wu=uwnd.2'
'define wv=vwnd.3'

****gradient***

'define r=6.37e6'
'define dtx = cdiff(pot,x)'
'define dty = cdiff(pot,y)'
'define dx = cdiff(lon,x)*r*cos(lat*3.1416/180)*3.1416/180'
'define dy = cdiff(lat,y)*r*3.1416/180'
'define ga=1/(mag(dtx/dx,dty/dy))'

******u,v winds************
'define dux=cdiff(wu,x)'
'define duy=cdiff(wu,y)'
'define dvx=cdiff(wv,x)'
'define dvy=cdiff(wv,y)'

'fr=-1*ga*((dtx/dx)*((dux/dx)*(dtx/dx)+(duy/dy)*(dty/dy))+(dty/dy)*((dvx/dx)*(dtx/dx)+(dvy/dy)*(dty/dy)))'
'm=m+fr'
tt=tt+1
endwhile
'm=m/15'

'd m'
;


改写成grd时只是在前后加了fwrite语句


'set gxout fwrite'
'set fwrite Z:\prog\grd\fro\try.grd'

'sdfopen z:\Guo\climate\pottem\pottmp.sig995.1991.nc'
'sdfopen Z:\Guo\climate\surfaceu\uwnd.sig995.1991.nc'
'sdfopen Z:\Guo\climate\surfacev\vwnd.sig995.1991.nc'

'set x 33 55'
'set y 41 59'

tt=225
'define m=0.0'
while (tt<=239)
'set t 'tt''
'define pot=pottmp.1'
'define wu=uwnd.2'
'define wv=vwnd.3'

****gradient***

'define r=6.37e6'
'define dtx = cdiff(pot,x)'
'define dty = cdiff(pot,y)'
'define dx = cdiff(lon,x)*r*cos(lat*3.1416/180)*3.1416/180'
'define dy = cdiff(lat,y)*r*3.1416/180'
'define ga=1/(mag(dtx/dx,dty/dy))'

******u,v winds************
'define dux=cdiff(wu,x)'
'define duy=cdiff(wu,y)'
'define dvx=cdiff(wv,x)'
'define dvy=cdiff(wv,y)'

*****frontogenesis******
'fr=-1*ga*((dtx/dx)*((dux/dx)*(dtx/dx)+(duy/dy)*(dty/dy))+(dty/dy)*((dvx/dx)*(dtx/dx)+(dvy/dy)*(dty/dy)))'
'm=m+fr'
tt=tt+1
endwhile

***avereage****
'm=m/15'

'd m'

'disable fwrite'
;



配合grd的ctl，我觉得问题可能在CTL上、直接拿别的文件的CTL改写的所以有些备注不太对哈，但是数据存放应该是对的。。

dset z:\prog\grd\fro\try.grd
title mean daily NMC reanalysis 1980-2013 June 1-30
undef -9.99e+33
xdef 23 linear 82.5 2.5
ydef 19 linear 10 2.5
zdef 1 levels 995
tdef 1 linear 00Z01Jun1980 1440mn
vars 1
m 0 99 mean frontogenesis
endvars

echo要读书

主楼里面我只截取了1991年的数据，可是实际我需要的是N年时间（1980-2013）不规律的数据，所以直接用nc画很麻烦，我就想着把N年的数据写到Grd里面再画图。

echo要读书

纠结在这老久了。。
请各位不吝赐教啊T T

river

你提取的不止一个时次吧，还有grads输出的资料默认缺测值都是-999000000.000000

278803532

你要确定经纬度的格点是不是从82.5和10开始的，我觉得你可以用'set lon'和‘set lat’这两个命令试试的，

echo要读书

river 发表于 2014-6-7 08:51
你提取的不止一个时次吧，还有grads输出的资料默认缺测值都是-999000000.000000

我提取的是比如t=222到t=333,求出每个时间的锋生函数之后又求了这段时间的平均锋生函数m,存放到文件里面的也是m，所以应该是一个t吧?

echo要读书

又试了下还是不行啊。。。是不是不能在转换nc-->grd的文件里面改写呢？。。但是不用这个改写没办法计算了。。。纠结。。请各位再帮我看看呢？附上使用的nc文件

echo要读书

问题和效果在主楼更新了一下，请各位帮忙看一下啊TAT @river @278803532


================06.09 补充一下问题========================
补充一下我想达到的效果以及方法。
我是想画出15年的特定时段（比如6-7月）平均锋生强度（根据文献，我将格点平均锋生强度叩定义为该格点上锋生函数为正的数值的累计值与锋生函数为正值事件的总频次的比值。），于是我在将nc转写grd文件过程中，先求出某年某时间段的平均锋生强度（但是在用if语句筛选的过程中貌似也出现了问题：if语句没有用？），再将15年的锋生函数平均。

由于我所需的时间年份跨度大（1991-2005），不能直接在gs文件里直接打开60多个文件，貌似我就必须得把所需的时间段的数据提到grd运算？同时又由于每年我所需的时间段不规则，所以我想直接在转换文件里面把平均锋生强度提出来写进grd文件。

故提取平均锋生强度的程序如下，主要出现的问题有两个：1.直接用nc进行运算的时候就没问题，改用grd的时候就会出现LZ那种问题；2.我用nc运算的时候，貌似if语句的筛选没有作用？画出来的图还是有负值的？


'set gxout fwrite'
'set fwrite Z:\prog\grd\fro\try.grd'

**********************************1991
'sdfopen z:\Guo\climate\temdaily\air.1991.nc'
'sdfopen z:\Guo\climate\rhumdaily\rhum.1991.nc'
'sdfopen Z:\Guo\climate\uwnd\uwnd.1991.nc'
'sdfopen Z:\Guo\climate\vwnd\vwnd.1991.nc'
'define all=0.0'

'set x 33 55'
'set y 41 59'
'set z 3'

tt=159
'define m=0.0'
'define day=0'
while (tt<=165)
'set t 'tt''
'define tem=air.1'
'define rh=rhum.2'
'define wu=uwnd.3'
'define wv=vwnd.4'

*****pot**********
'define prs=850'
'define es=(6.112*exp(17.67*(tem-273.15)/(tem-38)))'
'define q=rh*(0.62197*es/(prs-0.378*es))/100'
'define pot=tem*pow((1000/prs),(0.2854*(1.0-0.28*q)))'

****gradient***

'define r=6.37e6'
'define dtx = cdiff(pot,x)'
'define dty = cdiff(pot,y)'
'define dx = cdiff(lon,x)*r*cos(lat*3.1416/180)*3.1416/180'
'define dy = cdiff(lat,y)*r*3.1416/180'
'define ga=1/(mag(dtx/dx,dty/dy))'

******u,v winds************
'define dux=cdiff(wu,x)'
'define duy=cdiff(wu,y)'
'define dvx=cdiff(wv,x)'
'define dvy=cdiff(wv,y)'

'fr=-1*ga*((dtx/dx)*((dux/dx)*(dtx/dx)+(duy/dy)*(dty/dy))+(dty/dy)*((dvx/dx)*(dtx/dx)+(dvy/dy)*(dty/dy)))'

if(fr>=0)      ****此处if语句貌似没起作用？？？
'm=m+fr'
'day=day+1'
endif
tt=tt+1
endwhile

'd m/day'

;