LCOV - code coverage report
Current view: top level - common - mo_read_spatial_data.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 78 88 88.6 %
Date: 2024-02-27 14:40:25 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !> \file mo_read_spatial_data.f90
       2             : !> \brief \copybrief mo_read_spatial_data
       3             : !> \details \copydetails mo_read_spatial_data
       4             : 
       5             : !> \brief Reads spatial input data.
       6             : !> \details This module is to read spatial input data, e.g. dem, aspect, flow direction.
       7             : !! The module provides a subroutine for ASCII files.
       8             : !! (Subroutine for NetCDF files will come with release 5.1).
       9             : !! The data are read from the specified directory.
      10             : !> \authors Juliane Mai
      11             : !> \date Dec 2012
      12             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      13             : !! mHM is released under the LGPLv3+ license \license_note
      14             : !> \ingroup f_common
      15             : MODULE mo_read_spatial_data
      16             : 
      17             :   ! This module provides routines to read spatial data.
      18             : 
      19             :   ! Written  Juliane Mai, Jan 2013
      20             :   ! Modified
      21             : 
      22             :   USE mo_kind, ONLY : i4, dp
      23             :   USE mo_os, ONLY : check_path_isfile
      24             :   use mo_message, only: error_message
      25             : 
      26             :   IMPLICIT NONE
      27             : 
      28             :   PUBLIC :: read_header_ascii           ! Reads header of ASCII files
      29             :   PUBLIC :: read_spatial_data_ascii     ! Read ASCII  files
      30             :   ! PUBLIC :: read_spatial_data_nc      ! Read netCDF files -> will be implemented in release 5.1
      31             : 
      32             :   ! ------------------------------------------------------------------
      33             : 
      34             :   !    NAME
      35             :   !        read_spatial_data_ascii
      36             : 
      37             :   !    PURPOSE
      38             :   !>       \brief Reads spatial data files of ASCII format.
      39             : 
      40             :   !>       \details Reads spatial input data, e.g. dem, aspect, flow direction.
      41             : 
      42             :   !    HISTORY
      43             :   !>       \authors Juliane Mai
      44             : 
      45             :   !>       \date Jan 2013
      46             : 
      47             :   ! Modifications:
      48             :   ! Matthias Zink  Feb 2013 - , added interface and routine for datatype i4
      49             :   ! David Schaefer Mar 2015 - , removed double allocation of temporary data
      50             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
      51             : 
      52             :   INTERFACE  read_spatial_data_ascii
      53             :     MODULE PROCEDURE read_spatial_data_ascii_i4, read_spatial_data_ascii_dp
      54             :   END INTERFACE read_spatial_data_ascii
      55             : 
      56             :   ! ------------------------------------------------------------------
      57             : 
      58             :   PRIVATE
      59             : 
      60             :   ! ------------------------------------------------------------------
      61             : 
      62             : CONTAINS
      63             : 
      64             :   ! ------------------------------------------------------------------
      65             : 
      66             :   !    NAME
      67             :   !        read_spatial_data_ascii_dp
      68             : 
      69             :   !    PURPOSE
      70             :   !>       \brief TODO: add description
      71             : 
      72             :   !>       \details TODO: add description
      73             : 
      74             :   !    INTENT(IN)
      75             :   !>       \param[in] "character(len = *) :: filename" filename with location
      76             :   !>       \param[in] "integer(i4) :: fileunit"        unit for opening the file
      77             :   !>       \param[in] "integer(i4) :: header_nCols"    number of columns of data fields:
      78             :   !>       \param[in] "integer(i4) :: header_nRows"    number of rows of data fields:
      79             :   !>       \param[in] "real(dp) :: header_xllcorner"   header read in lower left corner
      80             :   !>       \param[in] "real(dp) :: header_yllcorner"   header read in lower left corner
      81             :   !>       \param[in] "real(dp) :: header_cellsize"    header read in cellsize
      82             : 
      83             :   !    INTENT(OUT)
      84             :   !>       \param[out] "real(dp), dimension(:, :) :: data" data
      85             :   !>       \param[out] "logical, dimension(:, :) :: mask"  mask
      86             : 
      87             :   !    HISTORY
      88             :   !>       \authors Robert Schweppe
      89             : 
      90             :   !>       \date Jun 2018
      91             : 
      92             :   ! Modifications:
      93             : 
      94          66 :   subroutine read_spatial_data_ascii_dp(filename, fileunit, header_ncols, header_nrows, header_xllcorner, &
      95             :                                        header_yllcorner, header_cellsize, data, mask)
      96             :     implicit none
      97             : 
      98             :     ! filename with location
      99             :     character(len = *), intent(in) :: filename
     100             : 
     101             :     ! unit for opening the file
     102             :     integer(i4), intent(in) :: fileunit
     103             : 
     104             :     ! number of rows of data fields:
     105             :     integer(i4), intent(in) :: header_nRows
     106             : 
     107             :     ! number of columns of data fields:
     108             :     integer(i4), intent(in) :: header_nCols
     109             : 
     110             :     ! header read in lower left corner
     111             :     real(dp), intent(in) :: header_xllcorner
     112             : 
     113             :     ! header read in lower left corner
     114             :     real(dp), intent(in) :: header_yllcorner
     115             : 
     116             :     ! header read in cellsize
     117             :     real(dp), intent(in) :: header_cellsize
     118             : 
     119             :     ! data
     120             :     real(dp), dimension(:, :), allocatable, intent(out) :: data
     121             : 
     122             :     ! mask
     123             :     logical, dimension(:, :), allocatable, intent(out) :: mask
     124             : 
     125             :     ! number of rows of data fields:
     126             :     integer(i4) :: file_nRows
     127             : 
     128             :     ! number of columns of data fields:
     129             :     integer(i4) :: file_nCols
     130             : 
     131             :     ! file read in lower left corner
     132          66 :     real(dp) :: file_xllcorner
     133             : 
     134             :     ! file read in lower left corner
     135          66 :     real(dp) :: file_yllcorner
     136             : 
     137             :     ! file read in cellsize
     138          66 :     real(dp) :: file_cellsize
     139             : 
     140             :     ! file read in nodata value
     141          66 :     real(dp) :: file_nodata
     142             : 
     143             :     integer(i4) :: i, j
     144             : 
     145             :     ! data
     146          66 :     real(dp), dimension(:, :), allocatable :: tmp_data
     147             : 
     148             :     ! mask
     149          66 :     logical, dimension(:, :), allocatable :: tmp_mask
     150             : 
     151             : 
     152             :     ! compare headers always with reference header (intent in)
     153             :     call read_header_ascii(filename, fileunit, &
     154             :             file_ncols, file_nrows, &
     155          66 :             file_xllcorner, file_yllcorner, file_cellsize, file_nodata)
     156          66 :     if ((file_ncols .ne. header_ncols)) &
     157           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: ncols')
     158          66 :     if ((file_nrows .ne. header_nrows)) &
     159           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: nrows')
     160          66 :     if ((abs(file_xllcorner - header_xllcorner) .gt. tiny(1.0_dp))) &
     161           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: xllcorner')
     162          66 :     if ((abs(file_yllcorner - header_yllcorner) .gt. tiny(1.0_dp))) &
     163           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: yllcorner')
     164          66 :     if ((abs(file_cellsize - header_cellsize)   .gt. tiny(1.0_dp))) &
     165           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: cellsize')
     166             : 
     167             :     ! allocation and initialization of matrices
     168         264 :     allocate(tmp_data(file_nrows, file_ncols))
     169     7829412 :     tmp_data = file_nodata
     170         264 :     allocate(tmp_mask(file_nrows, file_ncols))
     171     7829412 :     tmp_mask = .true.
     172             : 
     173             : 
     174             :     !checking whether the file exists
     175          66 :     call check_path_isfile(path = filename, raise=.true.)
     176             :     ! read in
     177             :     ! recl is only a rough estimate on bytes per line in the ascii
     178             :     ! default for nag: recl=1024(byte) which is not enough for 100s of columns
     179          66 :     open (unit = fileunit, file = filename, action = 'read', status = 'old', recl = 48 * file_ncols)
     180             :     ! (a) skip header
     181         462 :     do i = 1, 6
     182         462 :       read(fileunit, *)
     183             :     end do
     184             :     ! (b) read data
     185       27426 :     do i = 1, file_nrows
     186     7837986 :       read(fileunit, *) (tmp_data(i, j), j = 1, file_ncols)
     187             :     end do
     188          66 :     close(fileunit)
     189             : 
     190             :     ! set mask .false. if nodata value appeared
     191     7829346 :     where (abs(tmp_data - file_nodata) .lt. tiny(1.0_dp))
     192             :       tmp_mask = .false.
     193             :     end where
     194             : 
     195             :     ! transpose of data due to longitude-latitude ordering
     196         264 :     allocate(data(file_ncols, file_nrows))
     197     7838052 :     data = transpose(tmp_data)
     198          66 :     deallocate(tmp_data)
     199             : 
     200         264 :     allocate(mask(file_ncols, file_nrows))
     201     7838052 :     mask = transpose(tmp_mask)
     202          66 :     deallocate(tmp_mask)
     203             : 
     204         132 :   end subroutine read_spatial_data_ascii_dp
     205             : 
     206             :   !    NAME
     207             :   !        read_spatial_data_ascii_i4
     208             : 
     209             :   !    PURPOSE
     210             :   !>       \brief TODO: add description
     211             : 
     212             :   !>       \details TODO: add description
     213             : 
     214             :   !    INTENT(IN)
     215             :   !>       \param[in] "character(len = *) :: filename" filename with location
     216             :   !>       \param[in] "integer(i4) :: fileunit"        unit for opening the file
     217             :   !>       \param[in] "integer(i4) :: header_nCols"    number of columns of data fields:
     218             :   !>       \param[in] "integer(i4) :: header_nRows"    number of rows of data fields:
     219             :   !>       \param[in] "real(dp) :: header_xllcorner"   header read in lower left corner
     220             :   !>       \param[in] "real(dp) :: header_yllcorner"   header read in lower left corner
     221             :   !>       \param[in] "real(dp) :: header_cellsize"    header read in cellsize
     222             : 
     223             :   !    INTENT(OUT)
     224             :   !>       \param[out] "integer(i4), dimension(:, :) :: data" data
     225             :   !>       \param[out] "logical, dimension(:, :) :: mask"     mask
     226             : 
     227             :   !    HISTORY
     228             :   !>       \authors Robert Schweppe
     229             : 
     230             :   !>       \date Jun 2018
     231             : 
     232             :   ! Modifications:
     233             : 
     234         192 :   subroutine read_spatial_data_ascii_i4(filename, fileunit, header_ncols, header_nrows, header_xllcorner, &
     235             :                                        header_yllcorner, header_cellsize, data, mask)
     236             :     implicit none
     237             : 
     238             :     ! filename with location
     239             :     character(len = *), intent(in) :: filename
     240             : 
     241             :     ! unit for opening the file
     242             :     integer(i4), intent(in) :: fileunit
     243             : 
     244             :     ! number of rows of data fields:
     245             :     integer(i4), intent(in) :: header_nRows
     246             : 
     247             :     ! number of columns of data fields:
     248             :     integer(i4), intent(in) :: header_nCols
     249             : 
     250             :     ! header read in lower left corner
     251             :     real(dp), intent(in) :: header_xllcorner
     252             : 
     253             :     ! header read in lower left corner
     254             :     real(dp), intent(in) :: header_yllcorner
     255             : 
     256             :     ! header read in cellsize
     257             :     real(dp), intent(in) :: header_cellsize
     258             : 
     259             :     ! data
     260             :     integer(i4), dimension(:, :), allocatable, intent(out) :: data
     261             : 
     262             :     ! mask
     263             :     logical, dimension(:, :), allocatable, intent(out) :: mask
     264             : 
     265             :     ! number of rows of data fields:
     266             :     integer(i4) :: file_nRows
     267             : 
     268             :     ! number of columns of data fields:
     269             :     integer(i4) :: file_nCols
     270             : 
     271             :     ! file read in lower left corner
     272         192 :     real(dp) :: file_xllcorner
     273             : 
     274             :     ! file read in lower left corner
     275         192 :     real(dp) :: file_yllcorner
     276             : 
     277             :     ! file read in cellsize
     278         192 :     real(dp) :: file_cellsize
     279             : 
     280             :     ! file read in nodata value
     281         192 :     real(dp) :: file_nodata
     282             : 
     283             :     integer(i4) :: i, j
     284             : 
     285             :     ! data
     286         192 :     integer(i4), dimension(:, :), allocatable :: tmp_data
     287             : 
     288             :     ! mask
     289         192 :     logical, dimension(:, :), allocatable :: tmp_mask
     290             : 
     291             : 
     292             :     ! compare headers always with reference header (intent in)
     293             :     call read_header_ascii(filename, fileunit, &
     294             :             file_ncols, file_nrows, &
     295         192 :             file_xllcorner, file_yllcorner, file_cellsize, file_nodata)
     296         192 :     if ((file_ncols .ne. header_ncols)) &
     297           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: ncols')
     298         192 :     if ((file_nrows .ne. header_nrows)) &
     299           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: nrows')
     300         192 :     if ((abs(file_xllcorner - header_xllcorner) .gt. tiny(1.0_dp))) &
     301           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: xllcorner')
     302         192 :     if ((abs(file_yllcorner - header_yllcorner) .gt. tiny(1.0_dp))) &
     303           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: yllcorner')
     304         192 :     if ((abs(file_cellsize - header_cellsize)   .gt. tiny(1.0_dp))) &
     305           0 :              call error_message('read_spatial_data_ascii: header not matching with reference header: cellsize')
     306             : 
     307             :     ! allocation and initialization of matrices
     308         768 :     allocate(tmp_data(file_nrows, file_ncols))
     309    22673136 :     tmp_data = int(file_nodata, i4)
     310         768 :     allocate(tmp_mask(file_nrows, file_ncols))
     311    22673136 :     tmp_mask = .true.
     312             : 
     313             :     !checking whether the file exists
     314         192 :     call check_path_isfile(path = filename, raise=.true.)
     315             :     ! read in
     316             :     ! recl is only a rough estimate on bytes per line in the ascii
     317             :     ! default for nag: recl=1024(byte) which is not enough for 100s of columns
     318         192 :     open (unit = fileunit, file = filename, action = 'read', status = 'old', recl = 48 * file_ncols)
     319             :     ! (a) skip header
     320        1344 :     do i = 1, 6
     321        1344 :       read(fileunit, *)
     322             :     end do
     323             :     ! (b) read data
     324       79488 :     do i = 1, file_nrows
     325    22697856 :       read(fileunit, *) (tmp_data(i, j), j = 1, file_ncols)
     326             :     end do
     327         192 :     close(fileunit)
     328             : 
     329             :     ! set mask .false. if nodata value appeared
     330    22672944 :     where (tmp_data .EQ. int(file_nodata, i4))
     331             :       tmp_mask = .false.
     332             :     end where
     333             : 
     334             :     ! transpose of data due to longitude-latitude ordering
     335         768 :     allocate(data(file_ncols, file_nrows))
     336    22698048 :     data = transpose(tmp_data)
     337         192 :     deallocate(tmp_data)
     338             : 
     339         768 :     allocate(mask(file_ncols, file_nrows))
     340    22698048 :     mask = transpose(tmp_mask)
     341         192 :     deallocate(tmp_mask)
     342             : 
     343          66 :   end subroutine read_spatial_data_ascii_i4
     344             : 
     345             :   ! ------------------------------------------------------------------
     346             : 
     347             :   !    NAME
     348             :   !        read_header_ascii
     349             : 
     350             :   !    PURPOSE
     351             :   !>       \brief Reads header lines of ASCII files.
     352             : 
     353             :   !>       \details Reads header lines of ASCII files, e.g. dem, aspect, flow direction.
     354             : 
     355             :   !    INTENT(IN)
     356             :   !>       \param[in] "character(len = *) :: filename" Name of file and its location
     357             :   !>       \param[in] "integer(i4) :: fileunit"        File unit for open file
     358             : 
     359             :   !    INTENT(OUT)
     360             :   !>       \param[out] "integer(i4) :: header_nCols"  Reference number of columns
     361             :   !>       \param[out] "integer(i4) :: header_nRows"  Reference number of rows
     362             :   !>       \param[out] "real(dp) :: header_xllcorner" Reference lower left corner (x)
     363             :   !>       \param[out] "real(dp) :: header_yllcorner" Reference lower left corner (y)
     364             :   !>       \param[out] "real(dp) :: header_cellsize"  Reference cell size [m]
     365             :   !>       \param[out] "real(dp) :: header_nodata"    Reference nodata value
     366             : 
     367             :   !    HISTORY
     368             :   !>       \authors Juliane Mai
     369             : 
     370             :   !>       \date Jan 2013
     371             : 
     372             :   ! Modifications:
     373             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     374             : 
     375         306 :   subroutine read_header_ascii(filename, fileunit, header_ncols, header_nrows, header_xllcorner, header_yllcorner, &
     376             :                               header_cellsize, header_nodata)
     377         192 :     use mo_common_constants, only : nodata_dp
     378             :     implicit none
     379             : 
     380             :     ! Name of file and its location
     381             :     character(len = *), intent(in) :: filename
     382             : 
     383             :     ! File unit for open file
     384             :     integer(i4), intent(in) :: fileunit
     385             : 
     386             :     ! Reference number of rows
     387             :     integer(i4), intent(out) :: header_nRows
     388             : 
     389             :     ! Reference number of columns
     390             :     integer(i4), intent(out) :: header_nCols
     391             : 
     392             :     ! Reference lower left corner (x)
     393             :     real(dp), intent(out) :: header_xllcorner
     394             : 
     395             :     ! Reference lower left corner (y)
     396             :     real(dp), intent(out) :: header_yllcorner
     397             : 
     398             :     ! Reference cell size [m]
     399             :     real(dp), intent(out) :: header_cellsize
     400             : 
     401             :     ! Reference nodata value
     402             :     real(dp), intent(out) :: header_nodata
     403             : 
     404             :     character(5) :: dummy
     405             :     integer(i4) :: io
     406             : 
     407             :     !checking whether the file exists
     408         306 :     call check_path_isfile(path = filename, raise=.true.)
     409             :     ! reading header from a file
     410         306 :     open (unit = fileunit, file = filename, status = 'old')
     411         306 :     read (fileunit, *) dummy, header_nCols
     412         306 :     read (fileunit, *) dummy, header_nRows
     413         306 :     read (fileunit, *) dummy, header_xllcorner
     414         306 :     read (fileunit, *) dummy, header_yllcorner
     415         306 :     read (fileunit, *) dummy, header_cellsize
     416         306 :     read (fileunit, *, iostat=io) dummy, header_nodata
     417             :     ! EOF reached (nodata not present, use default value)
     418         306 :     if (io < 0) header_nodata = nodata_dp
     419         306 :     close(fileunit)
     420         306 :     dummy = dummy // ''   ! only to avoid warning
     421             : 
     422         306 :   end subroutine read_header_ascii
     423             : 
     424             : END MODULE mo_read_spatial_data

Generated by: LCOV version 1.16