黄乾

个人信息

Personal information

副教授     硕士生导师
性别:男
主要任职:大气物理系主任
在职信息:在岗
所在单位:大气物理学院
学科:大气科学

办公地点:气象楼1001

教师博客

[Grads画图]Fortran制作grads底图
发布时间:2021-07-22  点击次数:

第一步:准备地图边界数据,格式如下,其中6、2表示后续经纬度数据个数
6
119.908859252930, 29.101516723633
119.910301208496, 29.100149154663
119.909950256348, 29.098592758179
119.907882690430, 29.098100662231
119.907005310059, 29.096168518066
119.907432556152, 29.096067428589
2
119.903472900391, 29.094945907593
119.902725219727, 29.095489501953

第二步:做分段处理(生成grad地图,每段经纬度个数不能超过255个)
可用fortran程序处理
例:
c   处理段落超过255,分段      
  parameter(m=10000,mm=250)
real lat(m),lon(m)
integer aa,bb
integer k,k1
open(1,file="d:\map\map.txt")
open(2,file="d:\map\map.tmp")
do 10 j=1,m
read(1,*,end=100)aa
    do i=1,aa
    read(1,*)lon(i),lat(i)
    enddo
    if(aa<255) then
    k=1
       write(2,222)aa,"0"
       do l=1,aa
         write(2,111)lon(l),lat(l)
     enddo
    else
    k=aa/250
    k1=mod(aa,250)
      do l=1,k
      write(2,222)mm,"0"
       do n=1,250
      write(2,111)lon(250*(l-1)+n),lat(250*(l-1)+n)
   enddo
      enddo
       write(2,222)k1,"0"
       do ll=1,k1
       write(2,111)lon(250*k+ll),lat(250*k+ll)
    enddo
    endif
111   format(f11.7,",",f10.7)
222   format(i4,2x,a1)
10 enddo
100   write(*,*)"read over,"
end
第三步:生成地图文件
fortran程序处理
integer*1  ihead(3),ilon(4),ilat(4)
       integer lon,lat,npairs
       equivalence (lon,ilon),(lat,ilat)
       real alat,alon
       character(len=30) flname
       flname="map.tmp"
       open(1,file=flname,status='old')
       flname="map"
      open(3,file=flname, form='binary')
      ihead(1)=1
      ihead(2)=1
      nsegs=0
   99 read(1,*,end=100) npairs
      nsegs=nsegs+1
      if(npairs.gt.255) then
      print *,'Number of pairs in segment no.',nsegs,'exceeds 255 !'
      stop 'job aborted'
      endif
      ihead(3)=npairs
      irec=irec+1
      write(3)ihead
      write(6,*)(ihead(k),k=1,3),npairs
      do 10 i=1,npairs
       read(1,*) alon,alat
      lon=int(alon*10000.)+0.5
      lat=int((alat+90.)*10000.)+0.5
      write(3)(ilon(k),k=3,1,-1)
      write(3)(ilat(k),k=3,1,-1)
   10 continue
      go to 99
  100 print *,'Total number of segments processed= ',nsegs
      print *,irec
      close(1)
      close(2)
      close(3)
      stop ' Normal termination'
      end
第四步:绘制地图
将生成的map考到grads根目录下的dat文件夹中
gs文件中加入
'set mpdset map'
'draw map'
即可显示地图

关于去掉边界外的等值线处理,还需配合province-basemap.gs,该文件可在论坛里下
同时要制作map_out.txt文件
map_out.txt格式为:经度范围、纬度范围、draw、边界经纬度序列、经纬度范围左下角坐标、右下角坐标、右上角坐标、左上角坐标、左下角坐标。注意:边界经纬度序列必须是一连续闭合的。
map_out例:
119.0999985  120.9000015   28.3999996   29.7999992  draw  119.9088593  29.1015167  119.9103012  29.1001492………… 119.0999985  28.3999996  120.9000015  28.3999996  120.9000015  29.7999992  119.0999985  29.7999992  119.0999985  28.3999996