求画某市的气温分布图？（已解决，帮忙锁帖）

善人/jw

论坛里有站点插值画图的

mofangbao

gaofangjia 发表于 2011-11-23 16:07

                               
登录/注册后可看大图

试了好多方法，都不行，有没有人能帮我画了张。我附上我参考的一个文档。

哪里不行呢？你已经做到哪一步呢？

gaofangjia

mofangbao 发表于 2011-11-24 08:47

                               
登录/注册后可看大图

哪里不行呢？你已经做到哪一步呢？

用fortran转换数据，程序出错。所以后续的grads也没法做。
然后尝试使用surfer就行绘制，但是现在某市的bln底图还没有弄出来。

mofangbao

gaofangjia 发表于 2011-11-24 08:57

                               
登录/注册后可看大图

用fortran转换数据，程序出错。所以后续的grads也没法做。
然后尝试使用surfer就行绘制，但是现在某市的 ...

论坛有生成bln的方法，仔细找找哦

gaofangjia

感谢大家一路相伴，我已经找一个群里的同仁要了一份bln底图，并画出来了。再次感谢！

luoziwuhui

gaofangjia 发表于 2011-11-23 16:07

                               
登录/注册后可看大图

试了好多方法，都不行，有没有人能帮我画了张。我附上我参考的一个文档。

就是按照 那个年纪传的txt就是按照 那个年纪传的txt来弄，我试了一下。画出来的图如下：

另外 我想你的数据倒数第二个好像有问题 9948  18.23  159.59  25.7
  具体做法：

一 先把站点文件转化为二进制文件
program bintd
         parameter (st=3276) ! st means nums of station
         character stid(st)*8
     integer nflag,nlev
         real lat(st),lon(st),std(st),tim
         open(30,file='02070108.txt')
         open(40,file='td.dat',form='binary')
         do i=1,st
           read(30,'(a5,2f8.2,f7.1)')stid(i),lat(i),lon(i),std(i)
           write(*,*)stid(i),lat(i),lon(i),std(i)
         end do  
         close(30)
         tim=0.0
         nflag=1
         nlev=1
         do i=1,st
          write(40)stid(i),lat(i),lon(i),tim,nlev,nflag,std(i)
         end do
         nlev=0
         write(40)stid(3276),lat(3276),lon(3276),tim,nlev,nflag
         close(40)
         end
   
二 创建站点文件的描述文件  
  1单个时间变量 单个要素  td.ctl
dset td.dat
dtype station
stnmap td.map
undef 999.9
title daily td
tdef 1 linear Jan2002 1mo
vars 1
td 0 99 tdpoint data
endvars
四 生成相应的格点文件   wgrid.f90
program wgrid
        ! to get  grids
        real gtd(50,50)
        ! redefine the grids numbers when necessary
        open(10,file='/Models/czh/wgrid.dat',form='binary')
        do i=1,50
          do j=1,50
           gtd(i,j)=0
          end do
        end do
        write(10)((gtd(i,j),i=1,50),j=1,50)
        end
五 生成格点文件的描述文件 grid.ctl
dset F:\test\wgrid.dat
title gtd
undef -9.99E33
xdef 22 linear 109 0.04
ydef 22 linear 18 0.03
zdef 1 linear 1 1
tdef 1 linear Jan2010 1mo
vars 1
gtd 0 99 grid td
endvars
六 生成最后的gs文件 stat.gs
'reinit'
'open F:/test/grid.ctl'
'open F:/test/station0208.ctl'
'enable print F:/test/station0208.gmf'
'set mpdset cnworld'
'set lon 109 110'
'set lat 18.1 18.7'
'set t 1'
'define tds=oacres(gtd,temp.2)'
**'set cint 0.5'
'set gxout shaded'
'd tds'
'set gxout contour'
'd tds'
'print'
'printim F:/test/td.gif white'
'disable print'
'reinit'

gaofangjia

luoziwuhui 发表于 2011-11-28 11:18

                               
登录/注册后可看大图

就是按照 那个年纪传的txt就是按照 那个年纪传的txt来弄，我试了一下。画出来的图如下：

另外 我想你的 ...

确实，倒数第二个数据有问题。“9948 18.23 159.59 25.7”应改为“59948 18.23 109.59 25.7”，另外问一下，后面的程序应该怎么改？没有底图文件，果然画得没有surfer好看。

luoziwuhui

我觉得不是没有地图，是你的范围太小了，显示不出地图来，我用同样的方法差值画的降水就很好。
下面是画你的图的程序：
for 先把站点文件转化为二进制文件
   program bintd
         parameter (st=22) ! st means nums of station
         character stid(st)*8
          integer:: nflag,nlev
         real lat(st),lon(st),std(st),tim
         open(10,file='F:\test\station0208.txt')
         open(20,file='F:\test\station0208.dat',form='binary')
         do i=1,st
           read(10,'(a4,f7.2,f8.2,f6.1)')stid(i),lat(i),lon(i),std(i)
           write(*,*)stid(i),lat(i),lon(i),std(i)
           pause
         end do  
         close(10)
         tim=0.0
         nflag=1
         nlev=1
         do i=1,st
          write(20)stid(i),lat(i),lon(i),tim,nlev,nflag,std(i)
         end do
         nlev=0
         write(20)stid(22),lat(22),lon(22),tim,nlev,nflag
         close(20)
         end

生成相应的格点文件   wgrid.f90

     program wgrid
        ! to get  grids
        real gtd(22,22)
        ! redefine the grids numbers when necessary
        open(10,file='F:\test\wgrid.dat',form='binary')
        do i=1,22
          do j=1,22
           gtd(i,j)=0
          end do
        end do
        write(10)((gtd(i,j),i=1,22),j=1,22)
        end

ctl 自己写

gs
'reinit'
'open F:/test/grid.ctl'
'open F:/test/station0208.ctl'
'enable print F:/test/station0208.gmf'
'set mpdset cnworld'

'set lon 109 110'
'set lat 18.1 18.7'
'set t 1'
'define tds=oacres(gtd,temp.2)'
**'set cint 0.5'
'set gxout shaded'
'd tds'
'set gxout contour'
'd tds'
'print'
'printim F:/test/td.gif white'
'disable print'
'reinit'

gaofangjia

luoziwuhui 发表于 2011-12-3 10:38

                               
登录/注册后可看大图

我觉得不是没有地图，是你的范围太小了，显示不出地图来，我用同样的方法差值画的降水就很好。
下面是画你 ...

谢谢！收起来了放进程序库里，备用。

singapple

我也遇到同样的问题，学习一下~