close
프로필 배경
프로필 로고

해솔 | 해처럼 밝고 소나무처럼 바르게

Fortran · 2019. 8. 25. fullscreen 넓게보기

[Fortran] 포트란 지상 관측소를 기준으로 근접한 1/4/9 위성 픽셀에 대해 평균 수행

목차

 정보

  • 업무명    : 포트란 지상 관측소를 기준으로 근접한 1/4/9 위성 픽셀에 대해 평균 수행

  • 작성자    : 이상호

  • 작성일    : 2019-08-25

  • 설   명    :

  • 수정이력 :

 

 내용

[특징]

  • 위성 자료와 지상 관측소를 공간 일치를 위한 소프트웨어가 필요하며 이 프로그램은 이러한 목적을 이러한 목적을 달성하기 위해 고안된 소프트웨어

 

etc-image-0

 

[기능]

  • Terra/MODIS 자료에 대한 위경도 및 이진 자료 읽기

  • 해당 이진 자료를 이용하여 식생지수 (NDVI)로 변환

  • 지상 관측소의 위·경도를 기준으로 거리 계산 및 1/4/9 격자 평균 수행

 

[사용법]

  • 입력 자료 (Terra/MODIS 위·경도 및 이진 자료)를 동일 디렉터리에 저장

  • 소스 코드를 컴파일 (gfortran Average_Over_Satellite_Pixels_Relative_To_Ground_ Station.f90)

  • 소스 코드를 실행 (./a.out)

  • 각 관측소의 평균값을 확인

 

[사용 OS]

  • Linux

 

[사용 언어]

  • Fortran

 

 소스 코드

[전체]

      implicit none
      integer, parameter :: X_MAX = 2030 , Y_MAX = 1354 , XY_DATA = 2752681
      integer, parameter :: NaN = 0
!      real, parameter :: add_offset = 0.49 , scale_factor = 0.002 ! Emis
      real, parameter :: add_offset = 0.00 , scale_factor = 0.02    ! LST
      
      integer :: i, j, m, n, l, mm, ll, ii
      integer :: cc, cc1, cc2, cc3, cc4
      real :: read_lat, read_lon
      real :: MIN_dd, MIN_dd2
      real,allocatable,dimension( : , : ) :: NDVI,lat,lon,dd, dd2
      real,allocatable,dimension( : ) :: Point_lat, Point_lon
      allocate(lat(100000,100000), lon(100000,100000), NDVI(100000,100000), &
               dd(10000,10000), dd2(10000,10000), &
               Point_lat(10000), Point_lon(10000) )
      
      open(10,file='read_lat.dat',status='old',action='read')
      open(11,file='read_lon.dat',status='old',action='read')
      open(12,file='lat.dat')
      open(13,file='lon.dat')
      
      do i=1,XY_DATA
         read(10,'(f10.4)') read_lat
         read(11,'(f10.3)') read_lon

            if(read_lat /= 0.d0) then  ! 빈칸 = 0.0
            if(read_lat < 1.d0) read_lat=read_lat*(10.d0**4.d0) ! 정수는 1 이하로
            write(12,*) read_lat
            endif

            if(read_lon /= 0.d0) then
            if(read_lon < 1.d0) read_lon=read_lon*(10.d0**3.d0) 
            write(13,*) read_lon
            end if
      end do
      
      close (10) ; close (11) ; close (12) ; close (13)
      
      open(14,file='NDVI.dat',status='old',action='read')
      open(15,file='lat.dat')
      open(16,file='lon.dat')
      open(17,file='data.dat')
      open(18,file='Mean.dat')
       

      do i = 1, X_MAX
         read(14,*) NDVI(i,1:Y_MAX)
      enddo

      do i = 1, X_MAX
      do j = 1, Y_MAX
         read(15,*) lat(i,j)
         read(16,*) lon(i,j)
            if(NDVI(i,j) /= NaN) NDVI(i,j) = (NDVI(i,j)*scale_factor) + add_offset
         write(17,*) lat(i,j), lon(i,j), NDVI(i,j)
      end do
      end do
            
      Point_lat = reshape ( (/37.7515, 37.5714, 37.3376, 37.2723, 36.9704, 36.6392, &
                              36.3720, 36.0326, 36.0053, 35.8782, 35.8408, 35.5601, &
                              35.1729, 35.1047, 34.8169, 33.2938, 33.3868, 36.7796/), (/18/) )

      Point_lon = reshape ( (/128.8910, 126.9658, 127.9466, 126.9853, 127.9527, 127.4407, &
                              127.3721, 129.3796, 126.7614, 128.6522, 127.1172, 129.3203, &
                              126.8916, 129.0320, 126.3812, 126.1628, 126.8802, 127.1213/), (/18/) )

      ! 강릉(105), 서울(108), 원주(114), 수원(119), 충주(127), 청주(131)
      ! 대전(133), 포항(138), 군산(140), 대구(143), 전주(146), 울산(152)
      ! 광주(156), 부산(159), 목포(165), 고산(185), 성산(188), 천안(232)
	  do ii = 1, 18
            cc=0 ; cc1=0 ; cc2=0 ; cc3=0 ; cc4=0
       
            do i=1,X_MAX
            do j=1,Y_MAX
               dd(i,j)=sqrt(((Point_lon(ii)-lon(i,j))**2.d0)+((Point_lat(ii)-lat(i,j))**2.d0))
            end do 
            end do
      
            MIN_dd = minval( dd(1:X_MAX, 1:Y_MAX) )
                                                                     !  EX)   (305,244)   (305,245)   (305,246)
                                                                     !        (306,244)   (306,245)   (306,246)
                                                                     !        (307,244)   (307,245) 
            do i=1,X_MAX
            do j=1,Y_MAX
               if( MIN_dd == dd(i,j) ) then                          !  EX)   
                     cc = count( NDVI(i-1:i+1, j-1:j+1) /= NaN )     !                    (306,245)   
                     cc1 = count( NDVI(i-1:i, j-1:j) /= NaN )        !        
                     cc2 = count( NDVI(i-1:i, j:j+1) /= NaN )
                     cc3 = count( NDVI(i:i+1, j-1:j) /= NaN )
                     cc4 = count( NDVI(i:i+1, j:j+1) /= NaN )
      
                     do m=i-1,i+1
                     do l=j-1,j+1
                        if(mod(m+l,2)==mod(i+j,2) .and. l/=j) dd2(m,l)=sqrt(((Point_lon(ii)-lon(m,l))**2.d0)+((Point_lat(ii)-lat(m,l))**2.d0))
                     end do
                     end do
      
                     MIN_dd2 = minval( dd2(i-1:i+1, j-1:j+1), mask=dd2(i-1:i+1, j-1:j+1) /= 0.d0)
      
                                                                     !  EX)   (305,244)               (305,246)
                                                                     !                  cc1        cc2        
                                                                     ! 
                                                                     !                  cc3        cc4        
                                                                     !        (307,244)               (307,246)
                     do mm=i-1,i+1                               
                     do ll=j-1,j+1                               
                        if( MIN_dd2 == dd(mm,ll) ) then

!***************************************  TotalSUM(NDVI(i:i+1, j-1:j), mask=NDVI(i:i+1, j-1:j) /= NaN)/real(cc3), &
                  cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
      
      ! cc4
            if(mm>i-1 .and. ll>j-1) &
                  write(18,'(2i5,f10.4,3i5,f10.4,i5,f10.4)')  &
                  i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
                  mm, ll, cc4, SUM(NDVI(i:i+1, j:j+1), mask=NDVI(i:i+1, j:j+1) /= NaN)/real(cc4), &
                  cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)


!***************  Each ( fort.21, fort.22, fort.23, fort24 )  *************************************
      write(20+ii,*) 'ASOS Obs', Point_lat(ii), Point_lon(ii)
      write(20+ii,*) 'MODIS Obs', lat(i,j), lon(i,j)

      ! cc1
            if(mm<i+1 .and. ll<j+1) &
                  write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)')  &
                  i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
                  mm, ll, cc1, SUM(NDVI(i-1:i, j-1:j), mask=NDVI(i-1:i, j-1:j) /= NaN)/real(cc1), &
                  cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
      
      ! cc2
            if(mm<i+1 .and. ll>j-1) &
                  write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)')  &
                  i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
                  mm, ll, cc2, SUM(NDVI(i-1:i, j:j+1), mask=NDVI(i-1:i, j:j+1) /= NaN)/real(cc2), &
                  cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
      
      ! cc3
            if(mm>i-1 .and. ll<j+1) &
                  write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)')  &
                  i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
                  mm, ll, cc3, SUM(NDVI(i:i+1, j-1:j), mask=NDVI(i:i+1, j-1:j) /= NaN)/real(cc3), &
                  cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
      
      ! cc4
            if(mm>i-1 .and. ll>j-1) &
                  write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)')  &
                  i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
                  mm, ll, cc4, SUM(NDVI(i:i+1, j:j+1), mask=NDVI(i:i+1, j:j+1) /= NaN)/real(cc4), &
                  cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
!******************************************************************************************************      

                       end if
                                                                     
                     end do ! mm
                     end do ! jj
               end if
            end do  ! i
            end do  ! i
      enddo ! ii
 
      end

 

[GitHub GIST]

implicit none
integer, parameter :: X_MAX = 2030 , Y_MAX = 1354 , XY_DATA = 2752681
integer, parameter :: NaN = 0
! real, parameter :: add_offset = 0.49 , scale_factor = 0.002 ! Emis
real, parameter :: add_offset = 0.00 , scale_factor = 0.02 ! LST
integer :: i, j, m, n, l, mm, ll, ii
integer :: cc, cc1, cc2, cc3, cc4
real :: read_lat, read_lon
real :: MIN_dd, MIN_dd2
real,allocatable,dimension( : , : ) :: NDVI,lat,lon,dd, dd2
real,allocatable,dimension( : ) :: Point_lat, Point_lon
allocate(lat(100000,100000), lon(100000,100000), NDVI(100000,100000), &
dd(10000,10000), dd2(10000,10000), &
Point_lat(10000), Point_lon(10000) )
open(10,file='read_lat.dat',status='old',action='read')
open(11,file='read_lon.dat',status='old',action='read')
open(12,file='lat.dat')
open(13,file='lon.dat')
do i=1,XY_DATA
read(10,'(f10.4)') read_lat
read(11,'(f10.3)') read_lon
if(read_lat /= 0.d0) then ! 빈칸 = 0.0
if(read_lat < 1.d0) read_lat=read_lat*(10.d0**4.d0) ! 정수는 1 이하로
write(12,*) read_lat
endif
if(read_lon /= 0.d0) then
if(read_lon < 1.d0) read_lon=read_lon*(10.d0**3.d0)
write(13,*) read_lon
end if
end do
close (10) ; close (11) ; close (12) ; close (13)
open(14,file='NDVI.dat',status='old',action='read')
open(15,file='lat.dat')
open(16,file='lon.dat')
open(17,file='data.dat')
open(18,file='Mean.dat')
do i = 1, X_MAX
read(14,*) NDVI(i,1:Y_MAX)
enddo
do i = 1, X_MAX
do j = 1, Y_MAX
read(15,*) lat(i,j)
read(16,*) lon(i,j)
if(NDVI(i,j) /= NaN) NDVI(i,j) = (NDVI(i,j)*scale_factor) + add_offset
write(17,*) lat(i,j), lon(i,j), NDVI(i,j)
end do
end do
Point_lat = reshape ( (/37.7515, 37.5714, 37.3376, 37.2723, 36.9704, 36.6392, &
36.3720, 36.0326, 36.0053, 35.8782, 35.8408, 35.5601, &
35.1729, 35.1047, 34.8169, 33.2938, 33.3868, 36.7796/), (/18/) )
Point_lon = reshape ( (/128.8910, 126.9658, 127.9466, 126.9853, 127.9527, 127.4407, &
127.3721, 129.3796, 126.7614, 128.6522, 127.1172, 129.3203, &
126.8916, 129.0320, 126.3812, 126.1628, 126.8802, 127.1213/), (/18/) )
! 강릉(105), 서울(108), 원주(114), 수원(119), 충주(127), 청주(131)
! 대전(133), 포항(138), 군산(140), 대구(143), 전주(146), 울산(152)
! 광주(156), 부산(159), 목포(165), 고산(185), 성산(188), 천안(232)
do ii = 1, 18
cc=0 ; cc1=0 ; cc2=0 ; cc3=0 ; cc4=0
do i=1,X_MAX
do j=1,Y_MAX
dd(i,j)=sqrt(((Point_lon(ii)-lon(i,j))**2.d0)+((Point_lat(ii)-lat(i,j))**2.d0))
end do
end do
MIN_dd = minval( dd(1:X_MAX, 1:Y_MAX) )
! EX) (305,244) (305,245) (305,246)
! (306,244) (306,245) (306,246)
! (307,244) (307,245)
do i=1,X_MAX
do j=1,Y_MAX
if( MIN_dd == dd(i,j) ) then ! EX)
cc = count( NDVI(i-1:i+1, j-1:j+1) /= NaN ) ! (306,245)
cc1 = count( NDVI(i-1:i, j-1:j) /= NaN ) !
cc2 = count( NDVI(i-1:i, j:j+1) /= NaN )
cc3 = count( NDVI(i:i+1, j-1:j) /= NaN )
cc4 = count( NDVI(i:i+1, j:j+1) /= NaN )
do m=i-1,i+1
do l=j-1,j+1
if(mod(m+l,2)==mod(i+j,2) .and. l/=j) dd2(m,l)=sqrt(((Point_lon(ii)-lon(m,l))**2.d0)+((Point_lat(ii)-lat(m,l))**2.d0))
end do
end do
MIN_dd2 = minval( dd2(i-1:i+1, j-1:j+1), mask=dd2(i-1:i+1, j-1:j+1) /= 0.d0)
! EX) (305,244) (305,246)
! cc1 cc2
!
! cc3 cc4
! (307,244) (307,246)
do mm=i-1,i+1
do ll=j-1,j+1
if( MIN_dd2 == dd(mm,ll) ) then
!*************************************** TotalSUM(NDVI(i:i+1, j-1:j), mask=NDVI(i:i+1, j-1:j) /= NaN)/real(cc3), &
cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
! cc4
if(mm>i-1 .and. ll>j-1) &
write(18,'(2i5,f10.4,3i5,f10.4,i5,f10.4)') &
i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
mm, ll, cc4, SUM(NDVI(i:i+1, j:j+1), mask=NDVI(i:i+1, j:j+1) /= NaN)/real(cc4), &
cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
!*************** Each ( fort.21, fort.22, fort.23, fort24 ) *************************************
write(20+ii,*) 'ASOS Obs', Point_lat(ii), Point_lon(ii)
write(20+ii,*) 'MODIS Obs', lat(i,j), lon(i,j)
! cc1
if(mm<i+1 .and. ll<j+1) &
write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)') &
i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
mm, ll, cc1, SUM(NDVI(i-1:i, j-1:j), mask=NDVI(i-1:i, j-1:j) /= NaN)/real(cc1), &
cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
! cc2
if(mm<i+1 .and. ll>j-1) &
write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)') &
i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
mm, ll, cc2, SUM(NDVI(i-1:i, j:j+1), mask=NDVI(i-1:i, j:j+1) /= NaN)/real(cc2), &
cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
! cc3
if(mm>i-1 .and. ll<j+1) &
write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)') &
i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
mm, ll, cc3, SUM(NDVI(i:i+1, j-1:j), mask=NDVI(i:i+1, j-1:j) /= NaN)/real(cc3), &
cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
! cc4
if(mm>i-1 .and. ll>j-1) &
write(20+ii,'(2i5,f10.4,3i5,f10.4,i5,f10.4)') &
i, j, SUM(NDVI(i:i, j:j), mask=NDVI(i:i, j:j) /= NaN), &
mm, ll, cc4, SUM(NDVI(i:i+1, j:j+1), mask=NDVI(i:i+1, j:j+1) /= NaN)/real(cc4), &
cc, SUM(NDVI(i-1:i+1, j-1:j+1), mask=NDVI(i-1:i+1, j-1:j+1) /= NaN )/real(cc)
!******************************************************************************************************
end if
end do ! mm
end do ! jj
end if
end do ! i
end do ! i
enddo ! ii
end
view raw Average_Over_Satellite_Pixels_Relative_To_Ground_ Station.f90 hosted with ❤ by GitHub

 

 참고 문헌

[논문]

  • 없음

[보고서]

  • 없음

[URL]

  • 없음

 

 문의사항

[기상학/프로그래밍 언어]

  • sangho.lee.1990@gmail.com

[해양학/천문학/빅데이터]

  • saimang0804@gmail.com
저작자표시 비영리 변경금지 (새창열림)

'프로그래밍 언어 > Fortran' 카테고리의 다른 글

[Fortran] 포트란 퇴각검색을 이용한 스도쿠 풀이 알고리즘  (0) 2019.09.05
[Fortran, Gnuplot, ShellScript] 기상 자료를 이용한 통계 분석 및 가시화  (0) 2019.09.01
[Fortran] 포트란 끝말잇기 자가 학습 알고리즘  (2) 2019.07.28
[Fortran] 포트란 Arrey Function 예제  (0) 2019.07.26
[Fortran] 포트란 Himawari-8/AHI (2 km)을 CERES (20 km)로 공간 해상도 일치  (0) 2019.07.25
프로그래밍 언어/Fortran 관련 글
더 보기

티스토리툴바