swan中读取雷达拼图数据生成ctl和dat文件的Fortran源码

lnsyyzp

SWAN系统中的diamond 31类数据格式是二进制的三维雷达拼图数据格式,参考格式说明做了一个fortran程序,试用一下可以出图。数据说明完全用的是SWAN的系统说明书中的内容，特此作出版权说明。编译器是Compaq Visual Fortran 95。不做附件了，直接粘上。

program main
implicit none
integer i,j,k
character(len=100) :: filename_char,outname_char,ctlname_char,gsname_char,gmfname_char,pngname_char
character(len=12) ZonName; ! diamond 131 12个字节
character(len=38) DataName;!数据说明(例如 2008年5月19日雷达三维拼图)38个字节
character(len=8) Flag;  ! 文件标志，"swan"
character(len=8) Version;  ! 数据版本号，"1.0"
integer*2 year;!2008 两个字节
integer*2 month;!05  两个字节
integer*2 day;!19    两个字节
integer*2 hour,hourz;!14   两个字节
integer*2 minute;!31 两个字节
integer*2 interval ;  !两个字节        
integer*2 XNumGrids;!1300 两个字节
integer*2 YNumGrids;!800 两个字节
integer*2 ZNumGrids;!20  两个字节
integer*4 RadarCount; !拼图雷达数 四个字节
real*4 StartLon; !网格开始经度（左上角） 四个字节
real*4 StartLat; !网格开始纬度（左上角） 四个字节
real*4 CenterLon;!网格中心经度 四个字节
real*4 CenterLat;!网格中心纬度 四个字节
real*4 XReso; !经度方向分辨率 四个字节
real*4 YReso; !纬度方向分辨率 四个字节
real*4 ZhighGrids(40);!垂直方向的高度（单位km）数目根据ZnumGrids而得（最大40层） 160个字节。
character*16 RadarStationName(20);    !雷达站点名称, 20*16字节
real*4  RadarLongitude(20);      !雷达站点所在经度，单位：度， 4*20字节
real*4  RadarLatitude(20);       !雷达站点所在纬度，单位：度， 4*20字节
real*4  RadarAltitude(20);      !雷达所在海拔高度，单位：米， 4*20字节
integer*1    MosaicFlag(20);    !该雷达数据是否包含在本次拼图中，未包含:0，包含:1, 20字节
character*1    Reserved(172)
integer*1 dbz  ! 主要变量，变形后的雷达反射率因子，二进制整型，0-255，由原始观测值*2+66得到
character*6 tt
real origdbz  !原始观测值
write(*,*)'Input time for data: DDHHMM(Universial Time)'
read(*,'(a6)') tt
filename_char='e:\20110111\gradsout\Z_OTHE_RADAMOSAIC_201308'//tt//'.bin'
outname_char=filename_char(1:len_trim(filename_char)-4)//'.dat'
ctlname_char=filename_char(1:len_trim(filename_char)-4)//'.ctl'
gsname_char=filename_char(1:len_trim(filename_char)-4)//'.gs'
gmfname_char=filename_char(1:len_trim(filename_char)-4)//'.gmf'
pngname_char=filename_char(1:len_trim(filename_char)-4)//'.png'
open(200,file=filename_char,form='binary',access="SEQUENTIAL")
open(300,file=outname_char,form='binary',access="SEQUENTIAL")
read(200) ZonName
!write(*,*) ZonName
read(200) DataName
!write(*,*) DataName
read(200) Flag,Version,year,month,day,hour,minute,interval;
!write(*,*) Flag,Version,year,month,day,hour,minute,interval
read(200)XNumGrids,YNumGrids,ZNumGrids,RadarCount,StartLon,StartLat,CenterLon,CenterLat,XReso,YReso;
write(*,*) 'XNumGrids=',XNumGrids,'YNumGrids=',YNumGrids,'ZNumGrids=',ZNumGrids,RadarCount,StartLon,StartLat,CenterLon,CenterLat,XReso,YReso

do i=1,40
  read(200) ZhighGrids(i)
  write(*,*) ZhighGrids(i)
end do
do i=1,20
  read(200) RadarStationName(i)
end do
do i=1,20
  read(200) RadarLongitude(i)
end do
do i=1,20
  read(200) RadarLatitude(i)
end do
do i=1,20
  read(200) RadarAltitude(i)
end do
do i=1,20
  read(200) MosaicFlag(i)
  write(*,*) i,MosaicFlag(i)
end do
do i=1,172
  read(200) Reserved(i)
end do
do k=1,ZNumGrids
  do j=1,YNumGrids
    do i=1,XNumGrids
    read(200) dbz
if(dbz<0) then  !注意，读入的dbz如果是负值要做出转换，变成0-255
origdbz=(dbz+256-66.0)/2.0
else if(dbz==0) then
origdbz=0.0
else
origdbz=(dbz-66.0)/2.0
endif
write(300) origdbz
end do
  end do
end do
close(200)
close(300)
open(400,file=ctlname_char,status='unknown')
WRITE(400,'(A)')  'dset '//trim(outname_char)
WRITE(400,'(A)') 'title "Sample dbz data"'
WRITE(400,'(A)') 'options sequential yrev'
WRITE(400,'(A)') 'undef -99'
WRITE(400,'(A)') 'xdef 1200 linear 115.000000 0.010000'
WRITE(400,'(A)') 'ydef 1100 linear 35.000000 0.010000'
WRITE(400,'(A)') 'zdef 21 levels'
WRITE(400,'(10f8.1)') (ZhighGrids(i),i=1,ZNumGrids)
WRITE(400,'(A)')'tdef 1 linear 0Z11aug2011 1mo'
WRITE(400,'(A)')'vars 1'
WRITE(400,'(A)')'dbz 21  99 radar reflectivity factor'
WRITE(400,'(A)')'endvars'
close(400)
open(500,file=gsname_char,status='unknown')
WRITE(500,'(A)')"'"//'reinit'//"'"
WRITE(500,'(A)')"'"//'open '//trim(ctlname_char)//"'"
WRITE(500,'(A)')"'"//'enable print '//gmfname_char//"'"
WRITE(500,'(A)')"'"//'c'//"'"
WRITE(500,'(A)')"'"//'set display color white'//"'"
WRITE(500,'(A)')'*set mpdset hires cnworld'
WRITE(500,'(A)')"'"//'province-basemap liaon '//"'"
WRITE(500,'(A)')"'"//'cnbasemap'//"'"
WRITE(500,'(A)')"'"//'set lat 38.4 43.4'//"'"
WRITE(500,'(A)')"'"//'set lon 118.5 126'//"'"
WRITE(500,'(A)')"'"//'set z 6'//"'"
WRITE(500,'(A)')"'"//'set gxout shaded'//"'"
WRITE(500,'(A)')"'"//'set csmooth on'//"'"
WRITE(500,'(A)')'*让添色图更光滑*'
WRITE(500,'(A)')"'"//'set cmin 0'//"'"
WRITE(500,'(A)')"'"//'set grid off'//"'"
WRITE(500,'(A)')"'"//'d dbz'//"'"
WRITE(500,'(A)')"'"//'set ylpos 0 r'//"'"
WRITE(500,'(A)')"'"//'set cmin 0'//"'"
WRITE(500,'(A)')"'"//'d dbz'//"'"
WRITE(500,'(A)')"'"//'draw title '//filename_char(40:51)//"'"
WRITE(500,'(A)')"'"//'print'//"'"
WRITE(500,'(A)')"'"//'disable print'//"'"
WRITE(500,'(A)')"'"//'printim '//pngname_char//'  x600 y400 white'//"'"
WRITE(500,'(A)')"'"//"!"//'gxeps -i f:/20110111/gradsout/Z_OTHE_RADAMOSAIC_201101'//tt &
//'.gmf -o f:/20110111/gradsout/example'//tt//'.eps'//"'"
CLOSE(500)
stop
end program main

qixiangren34

放张图，让大家看看，我编译执行了，怎么出不来像样的图呢？

river

感谢分享，很棒

freekiller

非常感谢，顺求格式说明

lnsyyzp

7.1.        三维拼图数据格式
文件头格式，长度1024个字节
char ZonName[12];        // diamond 131 12个字节
char DataName[38];//数据说明(例如 2008年5月19日雷达三维拼图)38个字节
char        Flag[8];                // 文件标志，"swan"
char        Version[8];                // 数据版本号，"1.0"
unsigned short int year;//2008 两个字节
unsigned short int month;//05  两个字节
unsigned short int day;//19    两个字节
unsigned short int hour;//14   两个字节
unsigned short int minute;//31 两个字节
unsigned short int interval ;  //两个字节        
unsigned short int XNumGrids;//1300 两个字节
unsigned short int YNumGrids;//800 两个字节
unsigned short int ZNumGrids;//20  两个字节
int RadarCount; //拼图雷达数 四个字节
float StartLon; //网格开始经度（左上角） 四个字节
float StartLat; //网格开始纬度（左上角） 四个字节
float CenterLon;//网格中心经度 四个字节
float CenterLat;//网格中心纬度 四个字节
float XReso;        //经度方向分辨率 四个字节
float YReso;        //纬度方向分辨率 四个字节
float ZhighGrids[40];//垂直方向的高度（单位km）数目根据ZnumGrids而得（最大40层） 160个字节。
char RadarStationName[20][16];    //雷达站点名称,        20*16字节
float  RadarLongitude[20];      //雷达站点所在经度，单位：度， 4*20字节
float  RadarLatitude[20];       //雷达站点所在纬度，单位：度， 4*20字节
float  RadarAltitude[20];      //雷达所在海拔高度，单位：米， 4*20字节
unsigned char    MosaicFlag[20];    //该雷达数据是否包含在本次拼图中，未包含:0，包含:1, 20字节
char        Reserved[172];       
接下来是数据块，从底层到高层进行排列共ZnumGrids层。一个字节存储一个数据，值的范围0-255，2*dBZ+66等于该字节的值。每层的数据从起始点（左上角）开始，按维向（纬度y）减小写每行的经向（经度x增大）数据。

zoe_lover

{:eb502:}{:eb502:}{:eb502:}

yuanluochenai

感谢楼主

rossbyman

好多代码，好好看看

dqxiaoding

谢谢楼主分享

dqxiaoding

楼主你好，能分享一下TITAN的格式说明吗？

547168768

太棒了