โปรแกรมฟอร์แทรนการคำนวนระยะทางระหว่างพิกัดจุดภูมิศาสตร์สองจุด

โปรแกรมฟอร์แทรนการคำนวนระยะทางระหว่างพิกัดจุดภูมิศาสตร์สองจุด

2014.07.31

พอดีมีจุดพิกัดสองจุดที่ต้องการหาระยะทางเพื่อนำมาคำนวนผลในขั้นต่อไป มีสองวิธีที่เห็นเค้าใช้กันบ่อยๆคือ

1.     Spherical Law of Cosines

2.     Haversine Formula

หากข้อหลังมีความถูกต้องเยอะกว่า มีโค้ดให้เล่นเยอะแยะ http://rosettacode.org/wiki/Haversine_formula#Fortran ลองเล่นดูนะครับ พร้อมแนวคิด http://en.wikipedia.org/wiki/Haversine_formula

http://webtrough.wordpress.com/2011/09/19/latlon-distance-formula-in-excel-haversine-and-spherical-law-of-cosines/

program cal_dist_on_sphere

!-----------------------------------------------------------------------------------------------

! 2014.07.31

! Haversine codes :: http://rosettacode.org/wiki/Haversine_formula#Fortran

! Haversine formula :: http://www.ig.utexas.edu/outreach/googleearth/latlong.html

!------------------------------------------------------------------------------------------------

implicit none

real,parameter :: re=6372.8*1000.0    !- earth radius in m

real :: lat1,lat2                     !- lat lon in degree converted to radian

real :: lon1,lon2

real :: dlat,dlon                     !- difference in lat and lon in radian

real :: a,c                           !- Haversine parameter

real :: dist                          !- distance in meters

real, parameter :: pi = 4*atan(1.0)   !- exploit intrinsic atan to generate pi

real,parameter :: rd=pi/180.0         !- radian conversion degree to radian

lat1=15.0;lat2=16.0

lon1=100.0;lon2=100.0

dlat=(lat2-lat1)*rd

dlon=(lon2-lon1)*rd

lat1=lat1*rd

lat2=lat2*rd

!-Haversine formula

a=(sin(dlat/2))**2+cos(lat1)*cos(lat2)*(sin(dlon/2)**2)

c=2*atan2(sqrt(a),sqrt(1-a))

dist=re*c

print*,lat1,lon1,lat2,lon2,dist

endprogram cal_dist_on_sphere

ขั้นตอน-การเฉลี่ยแบบจำลองความสูงด้วย Kernel หรือ Moving window

2014.07.31

ขั้นตอน-การเฉลี่ยแบบจำลองความสูงด้วย Kernel หรือ Moving window

ผมต้องการจะเฉลี่ยค่าแบบจำลองความสูง GTOPO30 1ระวาง ด้วยพื้นที่โดยรอบของจุดที่พิจารณา ในที่นี้ใช้หลักการ  Moving window แล้วทดลองทำแต่ละระยะทางจากจุดความสูงที่กำลังพิจารณาที่ 10,25,50 เพื่อทำการเปรียบเทียบ และนำไปใช้ในงานอื่นต่อไป

1.                  ใช้ Fortran ในการเฉลี่ยตามโค้ดด้านล่าง

program avg_gtopo30

!----------------------------------------------------------

! 2014.07.30

! read gtopo 1 tile after swap to little endian

! apply average in horizontal length at considering grid

! write to binary and plot in grads

!----------------------------------------------------------

implicit none

character*150 :: infile_path

integer,parameter :: nrow=6000,ncol=4800                ! dem dimension

integer,parameter :: row=6000,col=4800

integer,parameter :: sizeofreal=4                       !-record size

integer,parameter :: sizeofint=2                        !-record size

integer*2,dimension(ncol,nrow) :: deml,demr                    !-dem left and right

real*4,dimension(col,row) :: demn                       !-merged dem

integer :: i,j,k

!-average dem in horizonal length

integer :: si,ei,sj,ej

integer :: m,n

real,parameter :: rgd=10.0                              !-prefered radius of grid [pixel]

real,dimension(ncol,nrow) :: demsm                      !-smooth dem output

real :: tdem    

real,parameter :: min_elv=-1.0E+8                       !- limit of minimum elevation                     

real,parameter :: max_elv= 1.0E+8                       !- limit of maximum elevation

real,parameter :: na_elv=-9999

infile_path='/home/mnattapon/RAID3storage4/research_icp/analysis_phb/beamblock_ans/dem_preparation/dem_raw_gtopo30/'

open(unit=10,file=trim(infile_path)//'E100N40.DEM2',form='unformatted',access='direct',recl=ncol*nrow*sizeofint)

print*,trim(infile_path)//'E100N40.DEM2'

read(10,rec=1)demr

close(10)

print*,'reading...'

do i=1,ncol

  do j=1,nrow

     m=ncol-i+1  !-col inverse

     n=nrow-j+1  !-row inverse

     demn(i,n)=demr(i,j)*1.0

  enddo

enddo

print*,'smooth...'

do i=1,ncol

  do j=1,nrow

     !-start and end pixel of defined radius

     si=i-rgd

     ei=i+rgd

     sj=j-rgd

     ej=j+rgd   

     !print*,i,j,si,ei,sj,ej

     !-unwanted area

     if ((i<=rgd) .or. ((ncol-i)<=rgd)) then

        demsm(i,j)=demn(i,j)

        cycle

     endif

     !-unwanted area

     if ((j<=rgd) .or. ((nrow-j)<=rgd)) then

        demsm(i,j)=demn(i,j)

        cycle

     endif

     !-loop for average dem

     tdem=0.0;k=0

     do m=si,ei

       do n=sj,ej

         !-check outlier in dem

         if ((demn(m,n)>max_elv).or.(demn(m,n)

             !print*,demn(m,n),'over!!'  

             cycle

         endif

         !print*,i,j,m,n,demn(m,n)

         tdem=tdem+demn(m,n)

         k=k+1

       enddo

     enddo 

     demsm(i,j)=tdem/k !-divided by total number of pixel

     !print*,demsm(i,j),demn(i,j),tdem,k

     !print*,'----------------------------------------------------------'

  enddo

enddo

open(unit=20,file='DEM_tile_E100N40_SM10.dat',form='unformatted',access='direct',recl=row*col*sizeofreal)

write(20,rec=1)demsm

close(20)

print*,'done...'

endprogram avg_gtopo30

2.                  สร้าง Control file ของ Grads ชื่อ ctl_read_dem1tile_sm10.ctl

dset DEM_tile_E100N40_SM10.dat

title GTOPO30 map 0.008333333333333 deg

undef -9999.0

**options little_endian template

**options little_endian

**options yrev little_endian

XDEF 4800  LINEAR   100     0.00833333333333

YDEF 6000  LINEAR  -10    0.00833333333333

zdef 1 levels 0 1

tdef 1 linear 22Z24aug2010 1hr

vars 1

dem     0 99 elevation [m]

Endvars

3.                  แสดงผลใน Grads

open ctl_read_dem1tile_sm10.ctl

set gxout shaded

d dem

cbarn

4.                  ตัวอย่างประกอบด้านล่าง