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
|