讨论：物理量随时间沿任意直线传播的剖面图画法

jolincai

背景：楼主想画这样的图，沿非平行于经纬线的任意直线（我已知道直线的起始点经纬度） 物理量随时间的变化图。如下所示，横坐标为经度和纬度，纵坐标为时间。楼主找了一份GS，可是对里面的有些算法实在是想不通，大家一起来看看解答一下！
GS如下：
'open e:\osci\Q1\jijie.ctl'
'set lon 60 140'
'set lat 0 50'    问题：这里的经纬度是我直线起始点的经纬度吗？
'set t 60'
*'define windspd=sqrt(u*u+v*v)'
'd ji'  
******以上是为了更好的确定经纬位置以及缩放比例********
say 'click a'  
'q bpos'
x0=subwrd(result,3)
y0=subwrd(result,4)   
*****以上是鼠标点击窗口确定剖线起始位置a（x0,y0）*******
say 'click b'
'q bpos'
x1=subwrd(result,3)
y1=subwrd(result,4)
******以上是鼠标点击窗口确定剖线终点位置b(x1,y1)******
'draw line 'x0' 'y0' 'x1' 'y1''   *在图形窗口中画出剖线*
'q xy2w 'x0' 'y0'' *将页面坐标（x0，y0）转换为对应的经纬坐标(lon0,lat0)*
lon0=subwrd(result,3)
lat0=subwrd(result,6)
'q xy2w 'x1' 'y1''*将页面坐标（x1，y1）转换为对应的经纬坐标（lon1，lat1）*
lon1=subwrd(result,3)
lat1=subwrd(result,6)
lx=(lon1-lon0)*112   *112是单位经度对应的实际距离（112km）*
ly=(lat1-lat0)*223    *223是单位纬度对应的实际距离（223km）*   重点问题：简单算一下 一个经度对应的实际距离也应该是周长除以360个经度=111KM啊？为什么一个纬度对应的实际距离成了223呢？不应该也是111吗？况且这个距离随纬度的不同而不同，高纬地区一个经度对应的实际距离只有赤道上的一半。这里直线的距离求法太粗糙了吧？
'd sqrt('lx'*'lx'+'ly'*'ly')'
l=subwrd(result,4)  *l是a，b两点的实际距离*
*decimal=l/55.8
decimal=l/50   *50是直线上相邻各点间的距离（水平分辨率）给出的值略大于模式水平分辨率为宜*  问题：为什么是除以50KM呢？和我资料的分辨率有关吗？2.5分辨率的也能除以50吗？
decimal_part=substr(decimal,3,20)
n=decimal-decimal_part  *这个确定decimal成为一个正数****
inclat=(lat1-lat0)/n   *inclat是相邻各等距点见的纬度增量*
inclon=(lon1-lon0)/n
****以下使计算a，b上个等分点所在经纬度******
i=1
while(i<=n)
lon.i=lon0+(i-1)*inclon
lat.i=lat0+(i-1)*inclat
i=i+1
endwhile
'set fwrite e:\osci\Q1\poumian.grd'
'set gxout fwrite'
nt=1
while(nt<=153)
i=1
'set t 'nt''
while(i<=n)
'set lon 'lon.i''
'set lat 'lat.i''
'd ji'
i=i+1
endwhile
nt=nt+1
endwhile
'disable fwrite'
fname = 'e:\osci\Q1\poumian.ctl'
rc=write(fname,'dset e:\osci\Q1\poumian.grd')
rc=write(fname,'undef -999.000')
rc=write(fname,'title Time Mask')
rc=write(fname,'xdef 'n' linear 'lon.1' 'lon.3' 'lon.5' 'lon.7' 'lon.9' 'lon.13' 'lon.16' 2.5')
rc=write(fname,'xdef 'n' linear 'lat.1' 'lat.3' 'lat.5' 'lat.7' 'lat.9' 'lat.13' 'lat.16' 2.5')
rc=write(fname,'ydef 1 linear 0 2.5')   问题：第一次接触XDEF既表示LAT又表示LON，既然把直线等分成了n份，那后面 'lon.1' 'lon.3' 'lon.5' 'lon.7' 'lon.9' 'lon.13' 'lon.16' 2.5'是什么意思，2.5表示什么？为什么不把所有等分点的lat lon都列出来？只列出第1 3 5 7 9 13 16呢？
rc=write(fname,'zdef 1 levels 1000')
rc=write(fname,'tdef 153  linear 01may2007 1dy')
rc=write(fname,'vars 1')
rc=write(fname,'pou 0 99 Time Mask')
rc=write(fname,'endvars')
rc=close(fname)
;
尽管有注释，楼主还是有疑问：问题为红色标注。

river

第一个问题，你设置的经纬度范围只要包括你要的直线的起始经纬度就可以。其实这个完全可以省略，直接给定起始点的经纬度就可以，屏幕上点击毕竟不会很精确。
第二个问题，建议使用grads的中央差分函数cdiff来分别计算dx 和 dy ，再利用地球半径来算出实际距离
第三个问题 同上
第四个问题，一看就是他写错了，应该是ydef ,而且他用的列举经纬度的办法，所以不能用linear，应该用levels，而且后面的2.5是不需要的。rc=write(fname,'ydef 1 linear 0 2.5') 这是多余的一句。而具体列举多少个经纬度的点，应该是和你提取过程中的那个 n 有关的吧。

jolincai

river 发表于 2017-6-7 16:53
第一个问题，你设置的经纬度范围只要包括你要的直线的起始经纬度就可以。其实这个完全可以省略，直接给定起 ...

问题2：单位经纬度对应的实际距离应该和分辨率没关系啊？不管分辨率多少，一个经纬度对应的距离都不变呀？这个脚本用勾股定理简单计算直线距离 可是我不懂112KM和223KM怎么来的。..
问题3：脚本在直线上每隔50KM取一个等分点，我就不知道我该取多少KM了。。。
然后写CTL时候，如果按你说的X Y T维都变化，那画图不又变成了横经度纵纬度的空间分布吗？脚本列举X方向，为什么不列举每一个，而是跳着列举呢？不应该从1到n点都列举出来么？
问题较多，大家一起讨论讨论吧~~~

jolincai

@路边摊1990 原创能来讨论一下吗？？

jolincai

@fatcat916 @edwardli 大神们来解疑答惑下~~~~·

jolincai

river 发表于 2017-6-7 16:53
第一个问题，你设置的经纬度范围只要包括你要的直线的起始经纬度就可以。其实这个完全可以省略，直接给定起 ...

另一个帖子里看到这个方法，简单明了，可是画出的图怪怪的，不知道脚本里set y 1是不是等于上一步中的lat值，如果不等于，那不相当于提取沿y 1的纬线上的数据了吗？
'reinit'
'sdfopen E:\uwnd.nc'
'set gxout fwrite '
'set fwrite E:\uwnd.dat'
'set grads off'
'set zlog on'
'set csmooth on'

'set t 1 20'
lon1 = 70.0
lon2 = 120.0
lat1 = 50.0
lat2 = 25.0
lon = lon1

while (lon <= lon2)
  lat = lat1 + (lat2-lat1)*(lon-lon1) / (lon2-lon1)
'set lat 'lat''
'set lon 'lon''
'a=u'
'set y 1'
'set lon 'lon''
'd a'
  lon = lon + 0.5

endwhile
'disable  fwrite '
'c'
'reinit '
;

river

jolincai 发表于 2017-6-8 17:12
另一个帖子里看到这个方法，简单明了，可是画出的图怪怪的，不知道脚本里set y 1是不是等于上一步中的lat ...

你看的那个帖子，别人就是要提取那个斜面上随时间变化的数据啊。我也提出了相同的疑问，但是楼主没有给我回复，我就忘记了

jolincai

help!!!!!!!!!

two_Drifters

请问这个叫什么图呀，想了解一下在matlab里画