分享一个自己写的【地理坐标转直角坐标小程序】，希望大神可以进来批评指正

little_five

最近在做台风路径突变的统计分析，需要先把经纬度转为平面坐标用的是中国台风1949-2015年的台风最佳路径集

看了好多资料和文献，发现真是隔行如隔山，什么高斯坐标、北京84，什么的完全不懂
后面搜到个简易方法，方法公式如下图

这步算出来的貌似没有问题
但是下一步要算台风路径的斜率和角度时，数据感觉很奇怪啊
求大神指导指导

最后附上我的程序

1. program main
   
2. implicit none
   
3.    
   
4.     character*5 :: liu
   
5.     character*4 :: wu,C
   
6.     integer :: i,j,k,line,ty_id
   
7.     character :: F,G
   
8.     character*20 :: name,h
   
9.     integer :: status = 0
   
10.     integer, allocatable :: ymdt(:),cata(:)
    
11.     real, allocatable :: lat(:),lon(:),pres(:),wind(:),x(:),y(:),w(:),N(:)
    
12.     character*8 :: year
    
13.     real*8, parameter :: pi=3.1415926
    
14.     real*8, parameter :: a=6378.137, b=6356.7523
    
15.     real*8 :: e2
    
16. 
    
17.     e2=(a**2-b**2)/a**2
    
18. 
    
19.     do k=1949,2015,1
    
20.       write(year,'(I4)')k
    
21.       open(unit = 10,file='h:/CH'//trim(adjustl(year))//'BST.txt')
    
22.       open(unit = 11,file='h:/CH'//trim(adjustl(year))//'BST.txt',status='replace',position='append')
    
23. 
    
24.         do while(.true.)
    
25.           read(10,*,iostat=status)liu,wu,line,ty_id,C,F,G,name,h
    
26.           if(status/=0)exit
    
27.           allocate(ymdt(line),cata(line),lat(line),lon(line),pres(line),wind(line),x(line),y(line),w(line),N(line))
    
28.           do i=1,line,1
    
29.             read(10,*,iostat=status)ymdt(i),cata(i),lat(i),lon(i),pres(i),wind(i)
    
30.             w(i)=dsqrt(1-((sin(lat(i)))**2)*e2)
    
31.             N(i)=a/w(i)
    
32.             x(i)=N(i)*cos(lat(i)/18000*pi)*cos(lon(i)/18000*pi)
    
33.             y(i)=N(i)*cos(lat(i)/18000*pi)*sin(lon(i)/18000*pi)
    
34.           end do
    
35.           write(11,"(A6,1X,I3,1X,I4,1X,A20)")liu,line,ty_id,name
    
36.           do j=1,line,1
    
37.             write(11,*)ymdt(j),cata(j),lat(j),lon(j),x(j),y(j)
    
38.           end do
    
39.        deallocate(ymdt,cata,lat,lon,pres,wind,x,y,w,N)
    
40.        end do
    
41.     end do
    
42. end

复制代码

ljchen1989

不是说西安80，北京54么？

little_five

ljchen1989 发表于 2018-1-28 16:08
不是说西安80，北京54么？

对对对，是这个的，但是我还是什么都不懂，然后看文献那个高斯-克吕格坐标转换，还有什么带之类的，也是看得一脸懵。。。。。。

little_five

ljchen1989 发表于 2018-1-28 16:08
不是说西安80，北京54么？

对对对，是这个的，但是我还是什么都不懂，然后看文献那个高斯-克吕格坐标转换，还有什么带之类的，也是看得一脸懵。。。。。。

931404656

已经不错啦  不过建议用proj4很强大

little_five

zhangqqqf 发表于 2018-1-29 07:58
觉得不错的样子，谢谢楼主

我的方法应该是错的，算出来的数据不对

little_five

931404656 发表于 2018-1-28 20:42
已经不错啦  不过建议用proj4很强大

我的方法应该是错的算出来的数据不对
然后搜了一天文献也没看出什么来
好郁闷{:cry:}

931404656

我用Fortran写过LCC转lonlat的

ljchen1989

little_five 发表于 2018-1-28 16:11
对对对，是这个的，但是我还是什么都不懂，然后看文献那个高斯-克吕格坐标转换，还有什么带之类的，也是 ...

据说要用3参数或5参数才能转换，然而这些参数属于国家机密
所以，从大地坐标到经纬度坐标，基本没有公式

little_five

931404656 发表于 2018-1-29 10:13
我用Fortran写过LCC转lonlat的

好厉害啊~
不过是我把问题想复杂了，现在不需要转了，哈哈