fortran提取任意区域内的点，可用来生成GrADS的mask文件

mofangbao

    今天看到@clare发帖说希望能提取平行四边形内的点，如果用在GrADS的maskout文件，这个用MeteoInfo实现实在很方便，如果有朋友希望自己通过fortran实现这个过程，那可以看看本帖子。
    最近在做的项目中也用到点和多边形的关系，于是就顺手转写为fortran的了，具体的算法大家可以查阅相关文献，这里我就不多说了。
    附件程序中我是利用这个函数输出了GrADS的maskout文件，并进行了测试，测试中我使用MeteoInfo提取中江苏省的边界（wmp文本格式），然后删掉了开始描述行，最后文件格式如图：


这里贴出生成GrADS的maskout文件的代码，关键子函数的代码请大家下载附件，就不收取贡献值了。

1. program GenMaskout
   
2. implicit none
   
3. integer,parameter::polys=1417,x=20,y=20                   !定义边界点和格点数量
   
4. real,parameter::interv=0.5,xs=115,ys=30                                                !定义格点开始经纬度和分辨率
   
5. integer,external::inRegion        
   
6. integer i,j
   
7. real pts(polys,2),pt(2),k                 !数组中1是横坐标  2是纵坐标
   
8. !读入你的边界点，polys是点数量，边界点一行一个点，经度在前，纬度在后，中间用逗号或者空格隔开
   
9. open(1,file='js.txt',status='old')         
   
10. do i=1,polys
    
11.   read(1,*)pts(i,1:2)
    
12. enddo
    
13. close(1)
    
14. open(1,file='js.grd',form='binary',status='replace')
    
15. pt(1)=xs
    
16. pt(2)=ys
    
17. do i=1,y
    
18.   do j=1,x
    
19.       k=inRegion(pts,pt,polys)   !这里调用判断点和多边形关系的函数
    
20.           write(1)k         
    
21.       pt(1)=pt(1)+interv
    
22.   enddo
    
23.   pt(2)=pt(2)+interv
    
24.   pt(1)=xs
    
25. enddo
    
26. close(1)
    
27. end
    
28. 
    

program GenMaskout<br /> implicit none<br /> integer,parameter::polys=1417,x=20,y=20                   !定义边界点和格点数量<br /> real,parameter::interv=0.5,xs=115,ys=30                                                !定义格点开始经纬度和分辨率<br /> integer,external::inRegion        <br /> integer i,j<br /> real pts(polys,2),pt(2),k                 !数组中1是横坐标  2是纵坐标<br /> !读入你的边界点，polys是点数量，边界点一行一个点，经度在前，纬度在后，中间用逗号或者空格隔开<br /> open(1,file='js.txt',status='old')          <br /> do i=1,polys<br />   read(1,*)pts(i,1:2)<br /> enddo<br /> close(1)<br /> open(1,file='js.grd',form='binary',status='replace')<br /> pt(1)=xs<br /> pt(2)=ys<br /> do i=1,y<br />   do j=1,x<br />       k=inRegion(pts,pt,polys)   !这里调用判断点和多边形关系的函数<br />           write(1)k          <br />       pt(1)=pt(1)+interv<br />   enddo<br />   pt(2)=pt(2)+interv<br />   pt(1)=xs<br /> enddo<br /> close(1)<br /> end<br /> <br />

最后在GrADS中测试一下吧！


这份代码的功能当然不仅仅是用来输出GrADS的maskout文件，其他用途大家可以尽情自由发挥，比如计算面雨量。。。

=======================
这里下载完整程序：
ptInRegion.f90 (1.99 KB, 下载次数: 133)

2013-11-14 16:47 上传

点击文件名下载附件

火_旦

lz fortran grads 玩的NB啊 膜拜~

polointer

膜拜大神！

风往北吹

玩转fortran真的接地气

风云翩跹

清风大大给力！

道狼52203

{:eb513:}{:eb513:}

kongfeng0824

清风出品，必属精品

kongfeng0824

清风出品，必属精品

心月培培

膜拜！！顶级！！MASTER！

LOSER

又学习了，谢谢！