自制grads地图

bekaku

第一步：准备地图边界数据，格式如下，其中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

如有问题，可邮件联系
bekaku@163.com




bekaku

自制的带乡镇边界底图，绘制的温度分布图效果

for_apply

bekaku 发表于 2013-6-5 15:11

                               
登录/注册后可看大图

自制的带乡镇边界底图，绘制的温度分布图效果

收下了，谢谢

loftydongshi

楼主辛苦了，给予诚挚的祝福！

liutaoz21

顶楼主

yl1992

很有用，谢谢

kongfeng0824

感谢！！！！！！！！！！！！

mofangbao

谢谢 这个是楼主自己写的还是网上的那个经典版本?这个程序应该是支持最基础的地图格式吧,想支持开关的那种楼主有没有试着尝试一下呢

qxtlyf

很有用，谢谢

bekaku

mofangbao 发表于 2013-6-6 08:58

                               
登录/注册后可看大图

谢谢 这个是楼主自己写的还是网上的那个经典版本?这个程序应该是支持最基础的地图格式吧,想支持开关的那种楼 ...

支持开关的能否介绍一下