Fortran으로 netCDF 파일 작성 (nf90_def_dim, def_var, put_att, put

Fortran

Fortran으로 netCDF 파일 작성 (nf90_def_dim, def_var, put_att, put_var)

이석사 중 2023. 8. 31. 15:57

728x90

이번 포스팅부터는 제목 앞에 작업이란 단어를 지우기로 했습니다

저는 나름대로 몇 가지 일을 했고 하나 하나가 제가 계획을 꾸려서 하는 작업이라고 생각했는데

논문을 써보면서 작업이라고 할 만한 것들이 아니라고 생각했습니다

앞으로는 제목을 간단하게 적고 중요했던 키포인트를 뒤쪽 괄호에 적는 식으로 제목을 바꿔보려고 합니다

오늘은 포트란으로 NetCDF 파일을 작성해보겠습니다

먼저 코드입니다

PROGRAM net2
  2         use netcdf
  3
  4         IMPLICIT NONE
  5
  6         INTEGER :: status, ncid
  7         INTEGER :: nx, ny, ii, jj, x, y, xx, yy
  8         REAL, DIMENSION(239, 239) :: field
  9         REAL, DIMENSION(239, 239) :: eqa
 10         INTEGER :: dimid_lon, dimid_lat, dimid_time
 11         INTEGER :: varid_lat, varid_lon, varid_field, varid_eqa, varid_time
 12         REAL :: pi, lat_array(239), lon_array(239)
 13         PARAMETER (pi = 4.0 * atan(1.0))
 14
 15         status = nf90_create('new.nc', NF90_NETCDF4, ncid)
 16         call check(status, 'open')
 17
 18         nx = 239
 19         ny = 239
 20
 21         status = nf90_def_dim(ncid, 'longitude', nx, dimid_lon)
 22         status = nf90_def_dim(ncid, 'latitude', ny, dimid_lat)
 23         status = nf90_def_dim(ncid, 'time', nf90_unlimited, dimid_time)
 24
 25         status = nf90_def_var(ncid, 'longitude', NF90_FLOAT, [dimid_lon], varid_lon)
 26         status = nf90_def_var(ncid, 'latitude', NF90_FLOAT, [dimid_lat], varid_lat)
 27         status = nf90_def_var(ncid, 'time', NF90_FLOAT, [dimid_time], varid_time)
 28         status = nf90_def_var(ncid, 'field', NF90_FLOAT, [dimid_lat, dimid_lon], varid_field)
 29         status = nf90_def_var(ncid, 'eqa', NF90_FLOAT, [dimid_lat, dimid_lon], varid_eqa)
 30
 31         status = nf90_put_att(ncid, NF90_GLOBAL, 'note', 'Test of Created NetCDF file')
 32         status = nf90_put_att(ncid, NF90_GLOBAL, 'ISIC', 15)
 33         status = nf90_put_att(ncid, varid_lon, 'units', 'degree_east')
 34         status = nf90_put_att(ncid, varid_lat, 'units', 'degree_north')
 35         status = nf90_put_att(ncid, varid_field, '_FillValue', -2e8)
 36         status = nf90_put_att(ncid, varid_eqa, '_FillValue', -2e8)
 37
 38         lat_array = [(jj * (360./ny), jj=0, ny-1)]
 39         lon_array = [((ii * (180./(nx-1)) - 90.), ii=0, nx-1)]
 40         do jj = 1, ny
 41                 do ii = 1, nx
 42                 field(ii, jj) = sin(lon_array(ii) * pi/180.) * cos(lat_array(jj) * pi/180.)
 43                 end do
 44         end do
 45
 46         do yy = 1, ny
 47                 do xx = 1, nx
 48                         eqa(xx, yy) = x**2 + y**2
 49                 end do
 50         end do
 51
 52         status = nf90_put_var(ncid, varid_lon, lon_array)
 53         status = nf90_put_var(ncid, varid_lat, lat_array)
 54         status = nf90_put_var(ncid, varid_field, field)
 55         status = nf90_put_var(ncid, varid_eqa, eqa)
 56
 57         status = nf90_close(ncid)
 58         call check(status, 'close')
 59
 60 contains
 61          SUBROUTINE check(status, loc)
 62                  INTEGER, INTENT(in) :: status
 63                  CHARACTER, INTENT(in) :: loc
 64                  if (status /= NF90_NOERR) then
 65                  WRITE(*, *) "ERROR at", loc
 66                  WRITE(*, *) NF90_STRERROR(status)
 67
 68                  end if
 69          END SUBROUTINE check
 70 END PROGRAM net2

중요한 것들을 하나씩 알아보겠습니다

  4         IMPLICIT NONE
  5
  6         INTEGER :: status, ncid
  7         INTEGER :: nx, ny, ii, jj, x, y, xx, yy
  8         REAL, DIMENSION(239, 239) :: field
  9         REAL, DIMENSION(239, 239) :: eqa
 10         INTEGER :: dimid_lon, dimid_lat, dimid_time
 11         INTEGER :: varid_lat, varid_lon, varid_field, varid_eqa, varid_time
 12         REAL :: pi, lat_array(239), lon_array(239)
 13         PARAMETER (pi = 4.0 * atan(1.0))

변수 선언부를 보면 굉장히 많은 변수들이 있습니다

더 중요한 부분은 dimid와 varid 변수들입니다

dimid는 차원에 대한 id이고 varid는 변수에 대한 id입니다

이 변수들은 integer로 선언했습니다

위도와 경도인 latitude와 longitude의 경우는 차원이 하나지만

위경도 격자 내에 입력되는 데이터는 시간축을 배제하면 위경도 2차원 격자를 가지기 때문에

field와 eqa는 DIMENSION을 이용해서 239 * 239 크기로 만들어줬습니다

 15         status = nf90_create('new.nc', NF90_NETCDF4, ncid)
 16         call check(status, 'open')
 17
 18         nx = 239
 19         ny = 239
 20
 21         status = nf90_def_dim(ncid, 'longitude', nx, dimid_lon)
 22         status = nf90_def_dim(ncid, 'latitude', ny, dimid_lat)
 23         status = nf90_def_dim(ncid, 'time', nf90_unlimited, dimid_time)
 24
 25         status = nf90_def_var(ncid, 'longitude', NF90_FLOAT, [dimid_lon], varid_lon)
 26         status = nf90_def_var(ncid, 'latitude', NF90_FLOAT, [dimid_lat], varid_lat)
 27         status = nf90_def_var(ncid, 'time', NF90_FLOAT, [dimid_time], varid_time)
 28         status = nf90_def_var(ncid, 'field', NF90_FLOAT, [dimid_lat, dimid_lon], varid_field)
 29         status = nf90_def_var(ncid, 'eqa', NF90_FLOAT, [dimid_lat, dimid_lon], varid_eqa)
 30
 31         status = nf90_put_att(ncid, NF90_GLOBAL, 'note', 'Test of Created NetCDF file')
 32         status = nf90_put_att(ncid, NF90_GLOBAL, 'ISIC', 15)
 33         status = nf90_put_att(ncid, varid_lon, 'units', 'degree_east')
 34         status = nf90_put_att(ncid, varid_lat, 'units', 'degree_north')
 35         status = nf90_put_att(ncid, varid_field, '_FillValue', -2e8)
 36         status = nf90_put_att(ncid, varid_eqa, '_FillValue', -2e8)

nf90_create로 먼저 netcdf 파일을 만들어줍니다

문법은 nf90_create(name, format, ncid) 형식입니다

다음은 nx, ny라는 변수에 격자의 개수를 지정해줬습니다

nf90_def_dim을 이용해서 먼저 차원을 선언해줍니다

제가 사용할 차원은 위도, 경도, 시간으로 총 3개입니다

위경도는 돌아가는 모델 내에서는 같기 때문에 잘 바꾸지 않고

시간의 경우는 1시간 간격, 3시간 간격처럼 변동이 있기 때문에 nf90_unlimited로 설정했습니다

다음은 nf90_def_var를 사용해서 변수를 선언해줍니다

저는 코드를 다른 웹사이트를 참고해서 할 때는 먼저 사이트의 코드를 그대로 입력시켜서 해보고

조금씩 수정하는 방식으로 하기 때문에 field라는 변수가 추가되어 있습니다

nf90_def_var의 문법은 nf90_def_var(ncid, name, format, dimension, varid) 형식입니다

nf90_def_ var는 사람으로 치면 개인정보에서 이름만 적어둔 것입니다

아직 나이, 생년월일 이런 정보들은 적히지 않은 것이죠

그래서 다음 nf90_put_att가 필요합니다

문법은 ndf90_put_att(ncid, varid, name, value) 형식입니다

이것들이 어떻게 추가되는지는 마지막에 나올 ncdump 결과를 보시면 더 이해가 빠르실 겁니다

 52 	    status = nf90_put_var(ncid, varid_lon, lon_array)
 53         status = nf90_put_var(ncid, varid_lat, lat_array)
 54         status = nf90_put_var(ncid, varid_field, field)
 55         status = nf90_put_var(ncid, varid_eqa, eqa)
 56
 57         status = nf90_close(ncid)
 58         call check(status, 'close')

앞에서 변수 이름과 속성들을 선언해줬지만 실제 데이터를 입력시키지 않았죠

nf90_put_var를 사용해서 데이터를 입력시켜줘야 합니다

n90_put_var(ncid, varid, data) 형식입니다

중간에 for loop를 사용한 부분이 임의의 값을 계산하는 부분입니다

마지막에 에러 판단을 위해 check 서브루틴을 넣어주면 코드는 끝입니다

마지막으로 ncdump 결과입니다

(base) lsh@DESKTOP-8N2HJ5V:~$ ncdump -hs new.nc
netcdf new {
dimensions:
        longitude = 239 ;
        latitude = 239 ;
        time = UNLIMITED ; // (0 currently)
variables:
        float longitude(longitude) ;
                longitude:units = "degree_east" ;
                longitude:_Storage = "contiguous" ;
                longitude:_Endianness = "little" ;
        float latitude(latitude) ;
                latitude:units = "degree_north" ;
                latitude:_Storage = "contiguous" ;
                latitude:_Endianness = "little" ;
        float time(time) ;
                time:_Storage = "chunked" ;
                time:_ChunkSizes = 1024 ;
                time:_Endianness = "little" ;
        float field(longitude, latitude) ;
                field:_FillValue = -2.e+08f ;
                field:_Storage = "contiguous" ;
                field:_Endianness = "little" ;
        float eqa(longitude, latitude) ;
                eqa:_FillValue = -2.e+08f ;
                eqa:_Storage = "contiguous" ;
                eqa:_Endianness = "little" ;

// global attributes:
                :note = "Test of Created NetCDF file" ;
                :ISIC = 15 ;
                :_NCProperties = "version=2,netcdf=4.9.0,hdf5=1.10.5" ;
                :_SuperblockVersion = 2 ;
                :_IsNetcdf4 = 1 ;
                :_Format = "netCDF-4" ;
}

longitude와 latitude를 보시면 nf90_put_att를 이용해서 units, 즉 단위를 입력시켰기 때문에

units = 'degree_east'라는 속성이 생겼습니다

추가로 global attributes라고 해서 전역 속성도 추가할 수 있습니다

1        status = nf90_put_att(ncid, NF90_GLOBAL, 'note', 'Test of Created NetCDF file')

위 처럼 varid가 들어갈 자리에 NF90_GLOBAL을 써주게 되면 추가가 가능합니다

728x90

'Fortran' 카테고리의 다른 글

Fortran으로 netCDF 파일 읽고 평균 구하기(gfortran 컴파일 에러) (0)	2023.08.24
작업7 : Fortran으로 파장 별 방출에너지 구하기 (0)	2022.10.26
작업6 : Fortran으로 합성 풍속 풍향 구하기2 (0)	2022.10.26
작업5 : Fortran으로 합성 풍향 풍속 구하기1 (0)	2022.10.26
작업4 : Fortran으로 시간방정식 계산하기 (0)	2022.10.26

현재글Fortran으로 netCDF 파일 작성 (nf90_def_dim, def_var, put_att, put_var)

대기과학지