Micaps第二类转GrADS的fortran源代码及其作图

mofangbao

   （谢谢射月楼主的提醒，该程序仅适用于类似EC的风场预报的格点数据）
    虽然我不常用这个，不过既然昨天看到群里有人需要，于是写出来供需要的参考一下。
代码很容易，新手可以自己试着转一下，代码中发现什么错误一定要指出来哦，嘿嘿
=====================================================
    先贴出fortran程序吧：

program main
integer,parameter::xnum=144,ynum=37
character*20 name
integer id,class
real lon,lat,hgt,pre,temp,dew_t,u(ynum,xnum),v(ynum,xnum)
!*****************变量说明*********************
! xnum:x格点数 ynum:y格点数                                        *
! name:存放读取以及生成新的文件的文件名                     *
! id,class,lon,lat,hgt,pre,temp.dew_t                                  *
! u(ynum,xnum),v(ynum,xnum)                                      *
! Micaps的风场的是以风速和风向来表示的，开始我以为 *
! 是像GrADS一样用u,v表示的，后来也没改变量名，       *
!介意的自己改下                                                             *
!**********************************************
!************读取需要转换的文件名（循环读取，每读取一次转换完再读第二个）*******************
open(1,file='name.txt',status='old',err=100)
do while(1)   
     read(1,*,end=200)name    !文件读完后跳出
     open(2,file=''//trim(name)//'',status='old',err=101)    !cvf为trim(name)即可/打开错误跳出
     name='b'//trim(name)
     print*,"正在转换-",trim(name),"-"
     open(10,file=''//trim(name)//'',status='replace',form='binary')    !建立新的二进制文件
     read(2,*)              !忽略文件头（两行）
     read(2,*)
     do i=1,ynum
          do j=1,xnum
               read(2,*)id,lon,lat,hgt,class,pre,temp,dew_t,u(i,j),v(i,j)   !循环读数，只有最后两个需要存储
          enddo
     enddo
201 close(2)
     do i=ynum,1,-1    !micaps文件存放顺序有点bt，和grads不一样，纬度从90到0的
          do j=1,xnum    !因此这里从最后一次写到第一组，便于grads处理
               write(10)u(i,j)    !写入风向
          enddo
     enddo
     do i=ynum,1,-1    !写入风速
          do j=1,xnum
               write(10)v(i,j)
          enddo
     enddo
     close(10)
enddo
goto 200
!到这里是因为读取出错了，给出提示并关闭文件、停止执行
100 print*, "文件name.txt读取错误！"
close(1)
stop
101 print*, "文件"//trim(name)//"读取错误！"
close(1)
close(2)
stop
200 pause "转换结束，按回车退出！"
close(1) !程序结束
end
=====================================================
    写这个fortran并不是一帆风顺的，主要是由于我对micaps何grads中风场的差别不太熟悉，后来弄清楚后，竟然发现，不说了，伤心啊...
    还是从百度百科找到了三角函数的相关公式才算法转换关系弄出来，结果发现原来是如此简单，转换可以在fortran中进行，也可以在grads中进行，我选择了在grads中进行，主要因为是不想再改这个程序了{:soso_e113:}
============================================================
由于本身就是格点数据，因此直接写了一个简单的测试ctl和gs文件，需要画图的朋友可能还要对图形进行一些美化，ctl文件如下：
dset E:\program\m2g\b11070720.000
title 850hpawind data
undef 9999
xdef 144 linear 0 2.5
ydef 37 linear 0 2.5
zdef 1 linear 1 1
tdef 1 linear jun2011 1dy
vars 2
u 0 99 u data
v 0 99 v data
endvars
=====ps:有兴趣的可以用批量描述试试
gs文件如下：
'reinit'
'open e:\program\m2g\m2g.ctl'
'set lon 60 160'
'set lat 0 70'
'set gxout barb'
*'d skip(-sin(u*3.14159/180)*v*2.5,2);skip(-cos(u*3.14159/180)*v*2.5,2);v'
*转换风速风向为U、V风场
'd -sin(u*3.14159/180)*v*2.5;-cos(u*3.14159/180)*v*2.5;v'
'draw title 11070720_000'
'printim e:\program\m2g\11070720_000.png white'
;
========================代码结束
好啦，给两张我测试的对比图，图形简陋，多多包涵啦
（这电脑装的是老古董-Micaps1，将就着点吧{:soso_e113:}）

11070720_024g

11070720_024m

==============================================

11070720_000g

11070720_000m

diamond2tograds.rar (863.17 KB, 下载次数: 420)

2011-7-8 20:05 上传

点击文件名下载附件

喜上眉梢

群主、楼主真厉害，非常非常感谢！！！

mofangbao

喜上眉梢 发表于 2011-7-8 20:09

                               
登录/注册后可看大图

群主、楼主真厉害，非常非常感谢！！！

不客气啦，互相学习

xuminge

为什么画出来的图中有好多圈圈，就像站点部分用圈表示，部分用点是吗？

mofangbao

xuminge 发表于 2011-7-9 13:26

                               
登录/注册后可看大图

为什么画出来的图中有好多圈圈，就像站点部分用圈表示，部分用点是吗？

说实话 那个圈圈我也不知道怎么会有的:L
等高手解释

xuminge

已经很感谢了，我正想要画风矢量的图了，谢谢！！！

majunlu

楼主真厉害

wangsansan

正在学习绘图，刚好学习学习，谢谢lz。

MeteoInfo

micaps第2类数据是高空全要素填图，应该是站点格式数据。micaps第11类格式是专门的格点矢量数据。不明白楼主为何用micaps第2类数据来表示格点矢量数据？

第十一类数据格式：格点矢量数据
文件头：
diamond  11  数据说明（字符串）  年  月  日  时次  时效  层次（均为整数）经度格距  纬度格距  起始经度  终止经度  起始纬度  终止纬度（均为浮点数）纬向格点数  经向格点数（均为整数）

注：此类数据主要用于画风场的流线。网格可以为经纬度网格，也可以为直角坐标网格。

当使用直角坐标网格时，文件头做如下改动：

1）把网格经度间隔和纬度间隔改为格点数据第一行最后一个点的经纬度。2）把起始经度和起始纬度改为格点数据第一行第一个点的经纬度。3）把终止经度和终止纬度改为格点数据最后一行最后一个点的经纬度。4）在第一行最后一点的经度上加一个数，指示在哪个底图投影下的直角坐标：LAMBERT投影加1000、MECATOR投影加2000、北半球投影加3000。

数据：

先放U分量，数据按先纬向后经向放（若为直角坐标网格数据，则先X方向，后Y方向），均为浮点数。所有格点的U分量放完后再放V分量，也是按先纬向后经向放。

mofangbao

MeteoInfo 发表于 2011-7-13 13:00

                               
登录/注册后可看大图

micaps第2类数据是高空全要素填图，应该是站点格式数据。micaps第11类格式是专门的格点矢量数据。不明白楼主 ...

谢谢MeteoInfo的指点
    从Micaps当初定义的数据格式来看，第二类要素的确是高空全要素填图，也当然是站点格式的数据，因为当初定义数据格式是这样的：

文件头：

diamond  2  数据说明（字符串）  年  月  日  时次  层次

总站点数（均为整数）

注：此类数据用于规范的高空填图

数据：

区站号（长整数）  经度  纬度  拔海高度（均为浮点数） 站点级别（整数）   高度  温度 温度露点差  风向  风速（均为浮点数）

=================
     所以我这样做是不通用的，当初写这个是因为群里有人就是想转换的EC的850hpa的风场预报，我也很奇怪为什么要把这个数据存为diamond2格式，因为这个EC的数据其实是2.5*2.5的格点数据，下面是数据片段
diamond 2 11070720_024时效EC_850hPa风场
    11     7     7    12   850  5328
   1     0.0    90.0    1    1 9999 9999 9999   74   11
   2     2.5    90.0    2    1 9999 9999 9999   76   11
   3     5.0    90.0    3    1 9999 9999 9999   79   11
   4     7.5    90.0    4    1 9999 9999 9999   81   11
   5    10.0    90.0    5    1 9999 9999 9999   84   11
   6    12.5    90.0    6    1 9999 9999 9999   86   11
   7    15.0    90.0    7    1 9999 9999 9999   89   11
   8    17.5    90.0    8    1 9999 9999 9999   91   11
   9    20.0    90.0    9    1 9999 9999 9999   94   11
  10    22.5    90.0   10    1 9999 9999 9999   96   11
  11    25.0    90.0   11    1 9999 9999 9999   99   11
  12    27.5    90.0   12    1 9999 9999 9999  101   11
  13    30.0    90.0   13    1 9999 9999 9999  104   11
  14    32.5    90.0   14    1 9999 9999 9999  106   11
  15    35.0    90.0   15    1 9999 9999 9999  109   11
  16    37.5    90.0   16    1 9999 9999 9999  111   11
  17    40.0    90.0   17    1 9999 9999 9999  114   11
================================
    可以看出，所谓的站号全部被赋值为类似格点ID（lon:144 lat:37），数据是从纬度最高处，经度最小处开始按照先经度后纬度的顺序开始排列，因此直接当做格点数据就可以写入grads进行画图，当然，如果是数据格式中定义的那种当然不能这么读取啦。
    所以这样读取只适合类似EC850风场预报的这种文件。
    我在其他地方也看过类似的讨论，说也可以按站点读入，然后再grads中进行插值，我还没有验证两者有什么区别。

    我上传了一个示例数据，射月楼主有兴趣的话可以试一下，Meteoinfo好像还没能直接读取diamond2格式，不知道我理解的对不对呢
11070720.024 (291.45 KB, 下载次数: 44)

2011-7-13 14:36 上传

点击文件名下载附件