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 : use mo_string_utils, only : num2str
26 :
27 : IMPLICIT NONE
28 : PUBLIC :: check_uniform_axis, get_header_info_from_nc, read_header_nc_or_ascii
29 : PUBLIC :: read_header_ascii ! Reads header of ASCII files
30 : PUBLIC :: read_spatial_data_nc_or_ascii ! Reads netcdf and ASCII files
31 : PUBLIC :: read_spatial_data_nc ! Reads netcdf files
32 : PUBLIC :: read_spatial_data_ascii ! Reads ASCII files
33 :
34 :
35 : ! ------------------------------------------------------------------
36 :
37 : ! NAME
38 : ! read_spatial_data_nc_or_ascii
39 :
40 : ! PURPOSE
41 : !> \brief Reads spatial data files of nc or ASCII format.
42 :
43 : !> \details Reads spatial input data, e.g. dem, aspect, flow direction.
44 :
45 : ! HISTORY
46 : !> \authors Simon Lüdke
47 :
48 : !> \date June 2025
49 :
50 : INTERFACE read_spatial_data_nc_or_ascii
51 : MODULE PROCEDURE read_spatial_data_nc_or_ascii_i4, read_spatial_data_nc_or_ascii_dp
52 : END INTERFACE read_spatial_data_nc_or_ascii
53 :
54 : ! ------------------------------------------------------------------
55 :
56 : ! NAME
57 : ! read_spatial_data_nc
58 :
59 : ! PURPOSE
60 : !> \brief Reads spatial data files of nc format.
61 :
62 : !> \details Reads spatial input data, e.g. dem, aspect, flow direction.
63 :
64 : ! HISTORY
65 : !> \authors Simon Lüdke
66 :
67 : !> \date June 2025
68 :
69 : INTERFACE read_spatial_data_nc
70 : MODULE PROCEDURE read_spatial_data_nc_i4, read_spatial_data_nc_dp
71 : END INTERFACE read_spatial_data_nc
72 :
73 : ! ------------------------------------------------------------------
74 :
75 : ! NAME
76 : ! read_spatial_data_ascii
77 :
78 : ! PURPOSE
79 : !> \brief Reads spatial data files of ASCII format.
80 :
81 : !> \details Reads spatial input data, e.g. dem, aspect, flow direction.
82 :
83 : ! HISTORY
84 : !> \authors Juliane Mai
85 :
86 : !> \date Jan 2013
87 :
88 : ! Modifications:
89 : ! Matthias Zink Feb 2013 - , added interface and routine for datatype i4
90 : ! David Schaefer Mar 2015 - , removed double allocation of temporary data
91 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
92 :
93 : INTERFACE read_spatial_data_ascii
94 : MODULE PROCEDURE read_spatial_data_ascii_i4, read_spatial_data_ascii_dp
95 : END INTERFACE read_spatial_data_ascii
96 :
97 : public :: check_header
98 :
99 :
100 :
101 : ! ------------------------------------------------------------------
102 :
103 : PRIVATE
104 :
105 : ! ------------------------------------------------------------------
106 :
107 : CONTAINS
108 :
109 : ! ------------------------------------------------------------------
110 :
111 : ! NAME
112 : ! read_spatial_data_ascii_dp
113 :
114 : ! PURPOSE
115 : !> \brief TODO: add description
116 :
117 : !> \details TODO: add description
118 :
119 : ! INTENT(IN)
120 : !> \param[in] "character(len = *) :: filename" filename with location
121 : !> \param[in] "integer(i4) :: fileunit" unit for opening the file
122 : !> \param[in] "integer(i4) :: header_nCols" number of columns of data fields:
123 : !> \param[in] "integer(i4) :: header_nRows" number of rows of data fields:
124 : !> \param[in] "real(dp) :: header_xllcorner" header read in lower left corner
125 : !> \param[in] "real(dp) :: header_yllcorner" header read in lower left corner
126 : !> \param[in] "real(dp) :: header_cellsize" header read in cellsize
127 :
128 : ! INTENT(OUT)
129 : !> \param[out] "real(dp), dimension(:, :) :: data" data
130 : !> \param[out] "logical, dimension(:, :) :: mask" mask
131 :
132 : ! HISTORY
133 : !> \authors Robert Schweppe
134 :
135 : !> \date Jun 2018
136 :
137 : ! Modifications:
138 :
139 66 : subroutine read_spatial_data_ascii_dp(filename, fileunit, header_ncols, header_nrows, header_xllcorner, &
140 : header_yllcorner, header_cellsize, data, mask)
141 : implicit none
142 :
143 : ! filename with location
144 : character(len = *), intent(in) :: filename
145 :
146 : ! unit for opening the file
147 : integer(i4), intent(in) :: fileunit
148 :
149 : ! number of rows of data fields:
150 : integer(i4), intent(in) :: header_nRows
151 :
152 : ! number of columns of data fields:
153 : integer(i4), intent(in) :: header_nCols
154 :
155 : ! header read in lower left corner
156 : real(dp), intent(in) :: header_xllcorner
157 :
158 : ! header read in lower left corner
159 : real(dp), intent(in) :: header_yllcorner
160 :
161 : ! header read in cellsize
162 : real(dp), intent(in) :: header_cellsize
163 :
164 : ! data
165 : real(dp), dimension(:, :), allocatable, intent(out) :: data
166 :
167 : ! mask
168 : logical, dimension(:, :), allocatable, intent(out) :: mask
169 :
170 : ! number of rows of data fields:
171 : integer(i4) :: file_nRows
172 :
173 : ! number of columns of data fields:
174 : integer(i4) :: file_nCols
175 :
176 : ! file read in lower left corner
177 66 : real(dp) :: file_xllcorner
178 :
179 : ! file read in lower left corner
180 66 : real(dp) :: file_yllcorner
181 :
182 : ! file read in cellsize
183 66 : real(dp) :: file_cellsize
184 :
185 : ! file read in nodata value
186 66 : real(dp) :: file_nodata
187 :
188 : integer(i4) :: i, j
189 :
190 : ! data
191 66 : real(dp), dimension(:, :), allocatable :: tmp_data
192 :
193 : ! mask
194 66 : logical, dimension(:, :), allocatable :: tmp_mask
195 :
196 : call read_header_ascii(filename, fileunit, &
197 : file_ncols, file_nrows, &
198 66 : file_xllcorner, file_yllcorner, file_cellsize, file_nodata)
199 : call check_header(file_ncols, file_nrows, file_xllcorner, file_yllcorner, file_cellsize, &
200 66 : header_ncols, header_nrows, header_xllcorner, header_yllcorner, header_cellsize)
201 : ! allocation and initialization of matrices
202 264 : allocate(tmp_data(file_nrows, file_ncols))
203 7829412 : tmp_data = file_nodata
204 264 : allocate(tmp_mask(file_nrows, file_ncols))
205 7829412 : tmp_mask = .true.
206 :
207 :
208 : !checking whether the file exists
209 66 : call check_path_isfile(path = filename, raise=.true.)
210 : ! read in
211 : ! recl is only a rough estimate on bytes per line in the ascii
212 : ! default for nag: recl=1024(byte) which is not enough for 100s of columns
213 66 : open (unit = fileunit, file = filename, action = 'read', status = 'old', recl = 48 * file_ncols)
214 : ! (a) skip header
215 462 : do i = 1, 6
216 462 : read(fileunit, *)
217 : end do
218 : ! (b) read data
219 27426 : do i = 1, file_nrows
220 7837986 : read(fileunit, *) (tmp_data(i, j), j = 1, file_ncols)
221 : end do
222 66 : close(fileunit)
223 :
224 : ! set mask .false. if nodata value appeared
225 7829346 : where (abs(tmp_data - file_nodata) .lt. tiny(1.0_dp))
226 : tmp_mask = .false.
227 : end where
228 :
229 : ! transpose of data due to longitude-latitude ordering
230 264 : allocate(data(file_ncols, file_nrows))
231 7838052 : data = transpose(tmp_data)
232 66 : deallocate(tmp_data)
233 :
234 264 : allocate(mask(file_ncols, file_nrows))
235 7838052 : mask = transpose(tmp_mask)
236 66 : deallocate(tmp_mask)
237 :
238 132 : end subroutine read_spatial_data_ascii_dp
239 :
240 : ! NAME
241 : ! read_spatial_data_ascii_i4
242 :
243 : ! PURPOSE
244 : !> \brief TODO: add description
245 :
246 : !> \details TODO: add description
247 :
248 : ! INTENT(IN)
249 : !> \param[in] "character(len = *) :: filename" filename with location
250 : !> \param[in] "integer(i4) :: fileunit" unit for opening the file
251 : !> \param[in] "integer(i4) :: header_nCols" number of columns of data fields:
252 : !> \param[in] "integer(i4) :: header_nRows" number of rows of data fields:
253 : !> \param[in] "real(dp) :: header_xllcorner" header read in lower left corner
254 : !> \param[in] "real(dp) :: header_yllcorner" header read in lower left corner
255 : !> \param[in] "real(dp) :: header_cellsize" header read in cellsize
256 :
257 : ! INTENT(OUT)
258 : !> \param[out] "integer(i4), dimension(:, :) :: data" data
259 : !> \param[out] "logical, dimension(:, :) :: mask" mask
260 :
261 : ! HISTORY
262 : !> \authors Robert Schweppe
263 :
264 : !> \date Jun 2018
265 :
266 : ! Modifications:
267 :
268 192 : subroutine read_spatial_data_ascii_i4(filename, fileunit, header_ncols, header_nrows, header_xllcorner, &
269 : header_yllcorner, header_cellsize, data, mask)
270 : implicit none
271 :
272 : ! filename with location
273 : character(len = *), intent(in) :: filename
274 :
275 : ! unit for opening the file
276 : integer(i4), intent(in) :: fileunit
277 :
278 : ! number of rows of data fields:
279 : integer(i4), intent(in) :: header_nRows
280 :
281 : ! number of columns of data fields:
282 : integer(i4), intent(in) :: header_nCols
283 :
284 : ! header read in lower left corner
285 : real(dp), intent(in) :: header_xllcorner
286 :
287 : ! header read in lower left corner
288 : real(dp), intent(in) :: header_yllcorner
289 :
290 : ! header read in cellsize
291 : real(dp), intent(in) :: header_cellsize
292 :
293 : ! data
294 : integer(i4), dimension(:, :), allocatable, intent(out) :: data
295 :
296 : ! mask
297 : logical, dimension(:, :), allocatable, intent(out) :: mask
298 :
299 : ! number of rows of data fields:
300 : integer(i4) :: file_nRows
301 :
302 : ! number of columns of data fields:
303 : integer(i4) :: file_nCols
304 :
305 : ! file read in lower left corner
306 192 : real(dp) :: file_xllcorner
307 :
308 : ! file read in lower left corner
309 192 : real(dp) :: file_yllcorner
310 :
311 : ! file read in cellsize
312 192 : real(dp) :: file_cellsize
313 :
314 : ! file read in nodata value
315 192 : real(dp) :: file_nodata
316 :
317 : integer(i4) :: i, j
318 :
319 : ! data
320 192 : integer(i4), dimension(:, :), allocatable :: tmp_data
321 :
322 : ! mask
323 192 : logical, dimension(:, :), allocatable :: tmp_mask
324 :
325 : call read_header_ascii(filename, fileunit, &
326 : file_ncols, file_nrows, &
327 192 : file_xllcorner, file_yllcorner, file_cellsize, file_nodata)
328 : call check_header(file_ncols, file_nrows, file_xllcorner, file_yllcorner, file_cellsize, &
329 192 : header_ncols, header_nrows, header_xllcorner, header_yllcorner, header_cellsize)
330 : ! allocation and initialization of matrices
331 768 : allocate(tmp_data(file_nrows, file_ncols))
332 22673136 : tmp_data = int(file_nodata, i4)
333 768 : allocate(tmp_mask(file_nrows, file_ncols))
334 22673136 : tmp_mask = .true.
335 :
336 : !checking whether the file exists
337 192 : call check_path_isfile(path = filename, raise=.true.)
338 : ! read in
339 : ! recl is only a rough estimate on bytes per line in the ascii
340 : ! default for nag: recl=1024(byte) which is not enough for 100s of columns
341 192 : open (unit = fileunit, file = filename, action = 'read', status = 'old', recl = 48 * file_ncols)
342 : ! (a) skip header
343 1344 : do i = 1, 6
344 1344 : read(fileunit, *)
345 : end do
346 : ! (b) read data
347 79488 : do i = 1, file_nrows
348 22697856 : read(fileunit, *) (tmp_data(i, j), j = 1, file_ncols)
349 : end do
350 192 : close(fileunit)
351 :
352 : ! set mask .false. if nodata value appeared
353 22672944 : where (tmp_data .EQ. int(file_nodata, i4))
354 : tmp_mask = .false.
355 : end where
356 :
357 : ! transpose of data due to longitude-latitude ordering
358 768 : allocate(data(file_ncols, file_nrows))
359 22698048 : data = transpose(tmp_data)
360 192 : deallocate(tmp_data)
361 :
362 768 : allocate(mask(file_ncols, file_nrows))
363 22698048 : mask = transpose(tmp_mask)
364 192 : deallocate(tmp_mask)
365 :
366 66 : end subroutine read_spatial_data_ascii_i4
367 :
368 : ! ------------------------------------------------------------------
369 :
370 : ! NAME
371 : ! read_header_ascii
372 :
373 : ! PURPOSE
374 : !> \brief Reads header lines of ASCII files.
375 :
376 : !> \details Reads header lines of ASCII files, e.g. dem, aspect, flow direction.
377 :
378 : ! INTENT(IN)
379 : !> \param[in] "character(len = *) :: filename" Name of file and its location
380 : !> \param[in] "integer(i4) :: fileunit" File unit for open file
381 :
382 : ! INTENT(OUT)
383 : !> \param[out] "integer(i4) :: header_nCols" Reference number of columns
384 : !> \param[out] "integer(i4) :: header_nRows" Reference number of rows
385 : !> \param[out] "real(dp) :: header_xllcorner" Reference lower left corner (x)
386 : !> \param[out] "real(dp) :: header_yllcorner" Reference lower left corner (y)
387 : !> \param[out] "real(dp) :: header_cellsize" Reference cell size [m]
388 : !> \param[out] "real(dp) :: header_nodata" Reference nodata value
389 :
390 : ! HISTORY
391 : !> \authors Juliane Mai
392 :
393 : !> \date Jan 2013
394 :
395 : ! Modifications:
396 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
397 :
398 306 : subroutine read_header_ascii(filename, fileunit, header_ncols, header_nrows, header_xllcorner, header_yllcorner, &
399 : header_cellsize, header_nodata)
400 192 : use mo_common_constants, only : nodata_dp
401 : implicit none
402 :
403 : ! Name of file and its location
404 : character(len = *), intent(in) :: filename
405 :
406 : ! File unit for open file
407 : integer(i4), intent(in) :: fileunit
408 :
409 : ! Reference number of rows
410 : integer(i4), intent(out) :: header_nRows
411 :
412 : ! Reference number of columns
413 : integer(i4), intent(out) :: header_nCols
414 :
415 : ! Reference lower left corner (x)
416 : real(dp), intent(out) :: header_xllcorner
417 :
418 : ! Reference lower left corner (y)
419 : real(dp), intent(out) :: header_yllcorner
420 :
421 : ! Reference cell size [m]
422 : real(dp), intent(out) :: header_cellsize
423 :
424 : ! Reference nodata value
425 : real(dp), intent(out) :: header_nodata
426 :
427 : character(5) :: dummy
428 : integer(i4) :: io
429 :
430 : !checking whether the file exists
431 306 : call check_path_isfile(path = filename, raise=.true.)
432 : ! reading header from a file
433 306 : open (unit = fileunit, file = filename, status = 'old')
434 306 : read (fileunit, *) dummy, header_nCols
435 306 : read (fileunit, *) dummy, header_nRows
436 306 : read (fileunit, *) dummy, header_xllcorner
437 306 : read (fileunit, *) dummy, header_yllcorner
438 306 : read (fileunit, *) dummy, header_cellsize
439 306 : read (fileunit, *, iostat=io) dummy, header_nodata
440 : ! EOF reached (nodata not present, use default value)
441 306 : if (io < 0) header_nodata = nodata_dp
442 306 : close(fileunit)
443 306 : dummy = dummy // '' ! only to avoid warning
444 :
445 306 : end subroutine read_header_ascii
446 :
447 284 : subroutine check_header(ncols, nrows, xllcorner, yllcorner, &
448 : cellsize, ref_ncols, ref_nrows, ref_xllcorner, ref_yllcorner, ref_cellsize, &
449 : tolerance, context)
450 :
451 306 : use mo_kind, only : dp, i4
452 : use mo_message, only : error_message
453 : use mo_string_utils, only : num2str
454 : use mo_common_constants, only: defaultTolerance_dp
455 :
456 : implicit none
457 :
458 : integer(i4), intent(in) :: ncols
459 : integer(i4), intent(in) :: nrows
460 : real(dp), intent(in) :: xllcorner
461 : real(dp), intent(in) :: yllcorner
462 : real(dp), intent(in) :: cellsize
463 : integer(i4), intent(in) :: ref_nrows
464 : integer(i4), intent(in) :: ref_ncols
465 : real(dp), intent(in) :: ref_xllcorner
466 : real(dp), intent(in) :: ref_yllcorner
467 : real(dp), intent(in) :: ref_cellsize
468 : real(dp), intent(in), optional :: tolerance
469 : character(*), intent(in), optional :: context
470 :
471 284 : real(dp) :: tol_local
472 284 : character(len=:), allocatable :: err_message
473 :
474 284 : if (present(tolerance)) then
475 0 : tol_local = tolerance
476 : else
477 : tol_local = defaultTolerance_dp
478 : end if
479 :
480 284 : err_message = ''
481 284 : if ((ref_ncols .ne. ncols)) then
482 : err_message = err_message // 'read_spatial_data_ascii: ncols mismatch file='// &
483 : trim(num2str(ncols))//' expected='//trim(num2str(ref_ncols))// &
484 0 : new_line('a')
485 : end if
486 284 : if ((nrows .ne. ref_nrows)) then
487 : err_message = err_message // 'read_spatial_data_ascii: nrows mismatch file='// &
488 : trim(num2str(nrows))//' expected='//trim(num2str(ref_nrows))// &
489 0 : new_line('a')
490 : end if
491 284 : if ((abs(xllcorner - ref_xllcorner) .gt. tol_local)) then
492 : err_message = err_message // 'read_spatial_data_ascii: xllcorner mismatch file='// &
493 : trim(num2str(xllcorner, '(G32.16)'))//' expected='// &
494 : trim(num2str(ref_xllcorner, '(G32.16)'))// &
495 0 : new_line('a')
496 : end if
497 284 : if ((abs(yllcorner - ref_yllcorner) .gt. tol_local)) then
498 : err_message = err_message // 'read_spatial_data_ascii: yllcorner mismatch file='// &
499 : trim(num2str(yllcorner, '(G32.16)'))//' expected='// &
500 : trim(num2str(ref_yllcorner, '(G32.16)'))// &
501 0 : new_line('a')
502 : end if
503 284 : if ((abs(cellsize - ref_cellsize) .gt. tol_local)) then
504 : err_message = err_message // 'read_spatial_data_ascii: cellsize mismatch file='// &
505 : trim(num2str(cellsize, '(G32.16)'))//' expected='// &
506 : trim(num2str(ref_cellsize, '(G32.16)'))// &
507 0 : new_line('a')
508 : end if
509 :
510 284 : if (len_trim(err_message) > 0) then
511 0 : if (present(context)) then
512 0 : err_message = trim(context) // new_line('a') // err_message
513 : end if
514 0 : call error_message(trim(err_message))
515 : end if
516 :
517 284 : end subroutine check_header
518 :
519 : !> \brief Read header information from either NetCDF or ASCII depending on file name/presence.
520 : !> \authors Simon Lüdke
521 : !> \date Feb 2025
522 22 : subroutine read_header_nc_or_ascii(filepath, fileunit, header_ncols, header_nrows, header_xllcorner, header_yllcorner, &
523 : header_cellsize, header_nodata, varName)
524 284 : use mo_os, only : path_root, path_isfile
525 : implicit none
526 : character(len=*), intent(in) :: filepath
527 : integer(i4), intent(in) :: fileunit
528 : integer(i4), intent(out) :: header_ncols, header_nrows
529 : real(dp), intent(out) :: header_xllcorner, header_yllcorner, header_cellsize, header_nodata
530 : character(len=*), intent(in), optional :: varName
531 :
532 22 : character(len=:), allocatable :: ncname
533 22 : character(len=:), allocatable :: trimmed
534 : integer(i4) :: len_path
535 : logical :: is_nc
536 :
537 22 : trimmed = trim(filepath)
538 22 : len_path = len_trim(trimmed)
539 22 : is_nc = .false.
540 22 : if (len_path >= 3) then
541 22 : if ((trimmed(len_path-2:len_path) == '.nc') .or. (trimmed(len_path-2:len_path) == '.NC')) is_nc = .true.
542 : end if
543 :
544 22 : ncname = trimmed
545 22 : if (.not. is_nc) ncname = path_root(trimmed) // '.nc'
546 :
547 44 : if (is_nc .or. path_isfile(ncname)) then
548 : call get_header_info_from_nc(ncname, fileunit, header_ncols, header_nrows, header_xllcorner, header_yllcorner, &
549 0 : header_cellsize, header_nodata, var=varName)
550 : else
551 22 : call read_header_ascii(trimmed, fileunit, header_ncols, header_nrows, header_xllcorner, header_yllcorner, &
552 : header_cellsize, header_nodata)
553 : end if
554 22 : end subroutine read_header_nc_or_ascii
555 :
556 : !> \brief check if given axis is a uniform axis.
557 : !> \authors Sebastian Müller
558 : !> \date Mar 2024
559 0 : subroutine check_uniform_axis(var, cellsize, origin, increasing, tol)
560 22 : use mo_netcdf, only : NcVariable
561 : use mo_utils, only: is_close
562 : implicit none
563 : type(NcVariable), intent(in) :: var !< NetCDF variable for corresponding axis
564 : real(dp), optional, intent(out) :: cellsize !< cellsize of the uniform axis
565 : real(dp), optional, intent(out) :: origin !< origin of the axis vertices
566 : logical, optional, intent(out) :: increasing !< whether the axis has increasing values
567 : real(dp), intent(in), optional :: tol !< tolerance for cell size comparisson (default: 1.e-7)
568 0 : real(dp), dimension(:), allocatable :: axis
569 0 : real(dp), dimension(:,:), allocatable :: bounds
570 0 : real(dp) :: diff, tol_
571 : integer(i4) :: i_ub, i_lb
572 : logical :: has_bnds
573 : type(NcVariable) :: bnds
574 : character(len=256) :: name
575 :
576 0 : call var%getData(axis)
577 0 : if (var%hasAttribute("bounds")) then
578 0 : has_bnds = .true.
579 0 : call var%getAttribute("bounds", name)
580 0 : bnds = var%parent%getVariable(trim(name))
581 0 : call bnds%getData(bounds)
582 : else
583 : has_bnds = .false.
584 : end if
585 : ! store var name for error messages
586 0 : name = var%getName()
587 :
588 0 : tol_ = 1.e-7_dp
589 0 : if ( present(tol) ) tol_ = tol
590 :
591 0 : if (size(axis) == 0_i4) &
592 0 : call error_message("check_uniform_axis: axis is empty: ", name)
593 :
594 0 : if (size(axis) > 1_i4) then
595 0 : diff = (axis(size(axis)) - axis(1)) / real(size(axis) - 1_i4, dp)
596 0 : if (.not.all(is_close(axis(2:size(axis))-axis(1:size(axis)-1), diff, rtol=0.0_dp, atol=tol_))) &
597 0 : call error_message("check_uniform_axis: given axis is not uniform: ", name)
598 : else
599 0 : if (.not. has_bnds) &
600 0 : call error_message("check_uniform_axis: can't check axis of size 1 when no bounds are given: ", name)
601 0 : diff = bounds(2,1) - bounds(1,1)
602 : end if
603 :
604 0 : if (has_bnds) then
605 : ! be forgiving if the bounds don't have the same monotonicity as the axis (cf-convetion is hard)
606 0 : i_lb = 1
607 0 : i_ub = 2
608 0 : if (size(bounds, dim=2)>1) then
609 0 : if (.not. is_close(bounds(2,1), bounds(1,2), rtol=0.0_dp, atol=tol_) &
610 0 : .and. is_close(bounds(1,1), bounds(2,2), rtol=0.0_dp, atol=tol_)) then
611 0 : print *, "check_uniform_axis: bounds actually have wrong monotonicity: ", name
612 0 : i_lb = 2
613 0 : i_ub = 1
614 : end if
615 : end if
616 0 : if (.not.all(is_close(bounds(i_ub,:)-bounds(i_lb,:), diff, rtol=0.0_dp, atol=tol_))) &
617 0 : call error_message("check_uniform_axis: given bounds are not uniform: ", name)
618 0 : if (.not.all(is_close(axis(:)-bounds(i_lb,:), 0.5_dp*diff, rtol=0.0_dp, atol=tol_))) &
619 0 : call error_message("check_uniform_axis: given bounds are not centered around axis points: ", name)
620 : end if
621 :
622 0 : if ( present(cellsize) ) cellsize = abs(diff)
623 0 : if ( present(origin) ) origin = minval(axis) - 0.5_dp * abs(diff)
624 0 : if ( present(increasing) ) increasing = diff > 0.0_dp
625 :
626 0 : end subroutine check_uniform_axis
627 :
628 : !> \brief check if given variable is a x-axis.
629 : !> \return `logical :: is_x_axis`
630 : !> \authors Sebastian Müller
631 : !> \date Mar 2024
632 0 : logical function is_x_axis(var)
633 0 : use mo_netcdf, only : NcVariable
634 : implicit none
635 : type(NcVariable), intent(in) :: var !< NetCDF variable to check
636 : character(len=256) :: tmp_str
637 :
638 0 : is_x_axis = .false.
639 0 : if (var%hasAttribute("standard_name")) then
640 0 : call var%getAttribute("standard_name", tmp_str)
641 0 : if (trim(tmp_str) == "projection_x_coordinate") is_x_axis = .true.
642 0 : else if (var%hasAttribute("axis")) then
643 0 : call var%getAttribute("axis", tmp_str)
644 0 : if (trim(tmp_str) == "X") is_x_axis = .true.
645 0 : else if (var%hasAttribute("_CoordinateAxisType")) then
646 0 : call var%getAttribute("_CoordinateAxisType", tmp_str)
647 0 : if (trim(tmp_str) == "GeoX") is_x_axis = .true.
648 : end if
649 0 : end function is_x_axis
650 :
651 : !> \brief check if given variable is a y-axis.
652 : !> \return `logical :: is_y_axis`
653 : !> \authors Sebastian Müller
654 : !> \date Mar 2024
655 0 : logical function is_y_axis(var)
656 0 : use mo_netcdf, only : NcVariable
657 : implicit none
658 : type(NcVariable), intent(in) :: var !< NetCDF variable to check
659 : character(len=256) :: tmp_str
660 :
661 0 : is_y_axis = .false.
662 0 : if (var%hasAttribute("standard_name")) then
663 0 : call var%getAttribute("standard_name", tmp_str)
664 0 : if (trim(tmp_str) == "projection_y_coordinate") is_y_axis = .true.
665 0 : else if (var%hasAttribute("axis")) then
666 0 : call var%getAttribute("axis", tmp_str)
667 0 : if (trim(tmp_str) == "Y") is_y_axis = .true.
668 0 : else if (var%hasAttribute("_CoordinateAxisType")) then
669 0 : call var%getAttribute("_CoordinateAxisType", tmp_str)
670 0 : if (trim(tmp_str) == "GeoY") is_y_axis = .true.
671 : end if
672 0 : end function is_y_axis
673 :
674 : !> \brief check if given variable is a lon coordinate.
675 : !> \return `logical :: is_lon_coord`
676 : !> \authors Sebastian Müller
677 : !> \date Mar 2024
678 0 : logical function is_lon_coord(var)
679 0 : use mo_netcdf, only : NcVariable
680 : implicit none
681 : type(NcVariable), intent(in) :: var !< NetCDF variable to check
682 : character(len=256) :: tmp_str
683 :
684 0 : is_lon_coord = .false.
685 0 : if (var%hasAttribute("standard_name")) then
686 0 : call var%getAttribute("standard_name", tmp_str)
687 0 : if (trim(tmp_str) == "longitude") is_lon_coord = .true.
688 0 : else if (var%hasAttribute("units")) then
689 0 : call var%getAttribute("units", tmp_str)
690 0 : if (trim(tmp_str) == "degreeE") is_lon_coord = .true.
691 0 : if (trim(tmp_str) == "degree_E") is_lon_coord = .true.
692 0 : if (trim(tmp_str) == "degree_east") is_lon_coord = .true.
693 0 : if (trim(tmp_str) == "degreesE") is_lon_coord = .true.
694 0 : if (trim(tmp_str) == "degrees_E") is_lon_coord = .true.
695 0 : if (trim(tmp_str) == "degrees_east") is_lon_coord = .true.
696 0 : else if (var%hasAttribute("_CoordinateAxisType")) then
697 0 : call var%getAttribute("_CoordinateAxisType", tmp_str)
698 0 : if (trim(tmp_str) == "Lon") is_lon_coord = .true.
699 0 : else if (var%hasAttribute("long_name")) then
700 0 : call var%getAttribute("long_name", tmp_str)
701 0 : if (trim(tmp_str) == "longitude") is_lon_coord = .true.
702 : end if
703 :
704 0 : end function is_lon_coord
705 :
706 : !> \brief check if given variable is a lat coordinate.
707 : !> \return `logical :: is_lat_coord`
708 : !> \authors Sebastian Müller
709 : !> \date Mar 2024
710 0 : logical function is_lat_coord(var)
711 0 : use mo_netcdf, only : NcVariable
712 : implicit none
713 : type(NcVariable), intent(in) :: var !< NetCDF variable to check
714 : character(len=256) :: tmp_str
715 :
716 0 : is_lat_coord = .false.
717 0 : if (var%hasAttribute("standard_name")) then
718 0 : call var%getAttribute("standard_name", tmp_str)
719 0 : if (trim(tmp_str) == "latitude") is_lat_coord = .true.
720 0 : else if (var%hasAttribute("units")) then
721 0 : call var%getAttribute("units", tmp_str)
722 0 : if (trim(tmp_str) == "degreeN") is_lat_coord = .true.
723 0 : if (trim(tmp_str) == "degree_N") is_lat_coord = .true.
724 0 : if (trim(tmp_str) == "degree_north") is_lat_coord = .true.
725 0 : if (trim(tmp_str) == "degreesN") is_lat_coord = .true.
726 0 : if (trim(tmp_str) == "degrees_N") is_lat_coord = .true.
727 0 : if (trim(tmp_str) == "degrees_north") is_lat_coord = .true.
728 0 : else if (var%hasAttribute("_CoordinateAxisType")) then
729 0 : call var%getAttribute("_CoordinateAxisType", tmp_str)
730 0 : if (trim(tmp_str) == "Lat") is_lat_coord = .true.
731 0 : else if (var%hasAttribute("long_name")) then
732 0 : call var%getAttribute("long_name", tmp_str)
733 0 : if (trim(tmp_str) == "latitude") is_lat_coord = .true.
734 : end if
735 0 : end function is_lat_coord
736 :
737 :
738 : !> \brief Read NetCDF header information (matching read_header_ascii signature)
739 : !> \details Reads grid metadata (ncols, nrows, corner, cellsize, nodata) from a NetCDF file.
740 : !> \authors Sebastian Müller, Simon Lüdke
741 : !> \date Mar 2024 / Feb 2025
742 0 : subroutine get_header_info_from_nc(ncname, fileunit, header_ncols, header_nrows, header_xllcorner, header_yllcorner, &
743 : header_cellsize, header_nodata, mask, var)
744 0 : use mo_netcdf, only : NcDataset, NcVariable, NcDimension
745 : use mo_utils, only : is_close
746 : use mo_string_utils, only : num2str
747 : implicit none
748 : character(*), intent(in) :: ncname !< NetCDF filename
749 : integer(i4), intent(in) :: fileunit !< Unused (kept for signature compatibility)
750 : integer(i4), intent(out) :: header_ncols !< number of columns
751 : integer(i4), intent(out) :: header_nrows !< number of rows
752 : real(dp), intent(out) :: header_xllcorner !< lower left x
753 : real(dp), intent(out) :: header_yllcorner !< lower left y
754 : real(dp), intent(out) :: header_cellsize !< cell size
755 : real(dp), intent(out) :: header_nodata !< nodata value
756 : logical, optional, allocatable, dimension(:,:), intent(out) :: mask !< mask
757 : character(*), intent(in), optional :: var !< optional variable name
758 :
759 : type(NcDataset) :: nc
760 : type(NcVariable) :: ncvar, xvar, yvar
761 0 : type(NcDimension), dimension(:), allocatable :: dims
762 : integer(i4) :: rnk, y_dir
763 : integer(i4) :: bottom_up, top_down
764 0 : real(dp) :: cs_x, cs_y, tol
765 : logical :: y_inc, x_sph, y_sph, x_cart, y_cart, flip_y
766 : integer(i4) :: keep_y
767 : character(len=256) :: varName
768 :
769 0 : keep_y = -1_i4
770 0 : y_dir = keep_y
771 :
772 : ! set defaults
773 0 : tol = 1.e-7_dp
774 0 : bottom_up = 1_i4
775 0 : top_down = 0_i4
776 :
777 0 : nc = NcDataset(ncname, "r")
778 :
779 0 : if (present(var)) then
780 0 : varName = trim(var)
781 : else
782 0 : varName = determine_data_var_name(nc, ncname)
783 : end if
784 :
785 0 : ncvar = nc%getVariable(varName)
786 0 : rnk = ncvar%getRank()
787 0 : if (rnk < 2) call error_message("grid % from_netcdf: given variable has too few dimensions: ", trim(nc%fname), ":", trim(varName))
788 :
789 0 : dims = ncvar%getDimensions()
790 0 : header_ncols = dims(1)%getLength()
791 0 : header_nrows = dims(2)%getLength()
792 0 : xvar = nc%getVariable(trim(dims(1)%getName()))
793 0 : yvar = nc%getVariable(trim(dims(2)%getName()))
794 :
795 : ! check if x/y axis are x/y/lon/lat by standard_name, units, axistype or long_name
796 0 : if (is_x_axis(yvar).or.is_lon_coord(yvar).or.is_y_axis(xvar).or.is_lat_coord(xvar)) &
797 0 : call error_message("grid % from_netcdf: variable seems to have wrong axis order (not y-x): ", trim(nc%fname), ":", trim(varName))
798 :
799 0 : x_cart = is_x_axis(xvar)
800 0 : y_cart = is_y_axis(yvar)
801 0 : x_sph = is_lon_coord(xvar)
802 0 : y_sph = is_lat_coord(yvar)
803 :
804 0 : if (.not.(x_cart.or.x_sph)) &
805 0 : call error_message("grid % from_netcdf: can't determine coordinate system from x-axis: ", trim(nc%fname), ":", trim(varName))
806 0 : if (.not.(y_cart.or.y_sph)) &
807 0 : call error_message("grid % from_netcdf: can't determine coordinate system from y-axis: ", trim(nc%fname), ":", trim(varName))
808 0 : if (.not.(x_sph.eqv.y_sph)) &
809 0 : call error_message("grid % from_netcdf: x and y axis seem to have different coordinate systems: ", trim(nc%fname), ":", trim(varName))
810 :
811 : ! check axis uniformity and monotonicity
812 0 : call check_uniform_axis(xvar, cellsize=cs_x, origin=header_xllcorner, tol=tol)
813 0 : call check_uniform_axis(yvar, cellsize=cs_y, origin=header_yllcorner, increasing=y_inc, tol=tol)
814 0 : if (y_dir == keep_y) then
815 0 : y_dir = top_down
816 0 : if (y_inc) y_dir = bottom_up
817 : end if
818 : ! check y_dir
819 0 : if (.not.any(y_dir==[bottom_up, top_down])) &
820 0 : call error_message("grid % from_netcdf: y-direction not valid: ", trim(num2str(y_dir)))
821 :
822 : ! warn about flipping if present axis is not in desired direction
823 0 : flip_y = y_inc.neqv.(y_dir==bottom_up)
824 0 : if (flip_y) then
825 0 : print *, "grid % from_netcdf: y axis direction is oposite to desired one (inefficient flipping). ", &
826 0 : "You could flip the file beforehand with: 'cdo invertlat <ifile> <ofile>'. ", trim(nc%fname), ":", trim(varName)
827 : end if
828 : ! check cellsize in x and y direction
829 0 : if (.not.is_close(cs_x, cs_y, rtol=0.0_dp, atol=tol)) &
830 0 : call error_message("grid % from_netcdf: x and y axis have different cell sizes: ", trim(nc%fname), ":", trim(varName))
831 0 : header_cellsize = cs_x
832 :
833 : ! Determine no-data value
834 0 : if (ncvar%hasAttribute("_FillValue")) then
835 0 : call ncvar%getAttribute("_FillValue", header_nodata)
836 0 : else if (ncvar%hasAttribute("missing_value")) then
837 0 : call ncvar%getAttribute("missing_value", header_nodata)
838 : else
839 0 : call error_message('***ERROR: read_nc_header: missing _FillValue or missing_value attribute')
840 : end if
841 :
842 0 : call nc%close()
843 0 : end subroutine get_header_info_from_nc
844 :
845 0 : function determine_data_var_name(nc, ncname) result(varName)
846 0 : use mo_netcdf, only : NcDataset, NcVariable
847 : use mo_message, only : error_message
848 : use mo_os, only : path_stem
849 : class(NcDataset), intent(in) :: nc
850 : character(len=*), intent(in) :: ncname
851 : character(len=256) :: varName
852 :
853 0 : type(NcVariable), dimension(:), allocatable :: vars
854 0 : character(len=256), dimension(:), allocatable :: data_var_names
855 : integer(i4) :: i, data_var_count, data_var_index
856 : character(len=256) :: stem
857 :
858 0 : varName = ''
859 0 : vars = nc%getVariables()
860 :
861 0 : data_var_count = 0
862 0 : do i = 1, size(vars)
863 0 : if (vars(i)%getNoDimensions() >= 2) data_var_count = data_var_count + 1
864 : end do
865 :
866 0 : if (data_var_count == 0) then
867 0 : call error_message('read_spatial_data_nc: could not determine variable name for file: ' // trim(ncname))
868 : end if
869 :
870 0 : allocate(data_var_names(data_var_count))
871 0 : data_var_index = 0
872 0 : do i = 1, size(vars)
873 0 : if (vars(i)%getNoDimensions() >= 2) then
874 0 : data_var_index = data_var_index + 1
875 0 : data_var_names(data_var_index) = trim(vars(i)%getName())
876 : end if
877 : end do
878 :
879 0 : if (data_var_count == 1) then
880 0 : varName = trim(data_var_names(1))
881 0 : return
882 : end if
883 :
884 0 : stem = trim(path_stem(ncname))
885 0 : do i = 1, data_var_count
886 0 : if ((trim(data_var_names(i)) == stem) .or. (trim(data_var_names(i)) == 'data')) then
887 0 : varName = trim(data_var_names(i))
888 0 : return
889 : end if
890 : end do
891 :
892 0 : call error_message('read_spatial_data_nc: variable name could not be determined for file: ' // trim(ncname))
893 0 : end function determine_data_var_name
894 :
895 :
896 : ! NAME
897 : ! read_spatial_data_nc_i4
898 :
899 : ! PURPOSE
900 : !> \brief Read file from filepath. If there is a nc file with the correct name it is prefered over the asci file.
901 :
902 : !> \details TODO: add description
903 :
904 : ! HISTORY
905 : !> \authors Simon Lüdke
906 :
907 : !> \date June 2025
908 0 : subroutine read_spatial_data_nc_i4(ncname, data, maskout, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value, varName)
909 0 : use mo_netcdf, only : NcDataset, NcVariable
910 : use mo_message, only : error_message
911 : implicit none
912 :
913 : character(len=*), intent(in) :: ncname
914 : character(len=*), intent(in), optional :: varName
915 : integer(i4), dimension(:, :), allocatable, intent(out) :: data
916 : logical, optional, allocatable, dimension(:, :), intent(out) :: maskout
917 : integer(i4), intent(out) :: ncols, nrows
918 : real(dp), intent(out) :: xllcorner, yllcorner, cellsize
919 : real(dp), intent(out) :: nodata_value
920 :
921 : type(NcDataset) :: nc
922 : type(NcVariable) :: var
923 0 : integer(i4), allocatable :: var_shape(:), start(:), cnt(:)
924 : character(len=256) :: varName_local
925 :
926 : ! Open NetCDF dataset
927 0 : nc = NcDataset(ncname, "r")
928 :
929 0 : if (present(varName)) then
930 0 : varName_local = trim(varName)
931 : else
932 0 : varName_local = determine_data_var_name(nc, ncname)
933 : end if
934 :
935 : ! Extract header info
936 0 : call get_header_info_from_nc(ncname, 0_i4, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value, var=varName_local)
937 :
938 : ! Retrieve variable
939 0 : var = nc%getVariable(trim(varName_local))
940 :
941 : ! Determine shape and prepare read extents
942 0 : var_shape = var%getShape()
943 0 : allocate(start(size(var_shape)), source=1_i4)
944 0 : allocate(cnt(size(var_shape)), source=1_i4)
945 0 : cnt(1:2) = var_shape(1:2)
946 :
947 : ! Determine no-data value
948 0 : if (var%hasAttribute("_FillValue")) then
949 0 : call var%getAttribute("_FillValue", nodata_value)
950 0 : else if (var%hasAttribute("missing_value")) then
951 0 : call var%getAttribute("missing_value", nodata_value)
952 : else
953 0 : call error_message('***ERROR: read_nc_i4_data: missing _FillValue or missing_value attribute')
954 : end if
955 :
956 : ! Allocate and read data
957 0 : allocate(data(var_shape(1), var_shape(2)))
958 0 : call var%getData(data, start=start, cnt=cnt)
959 :
960 0 : if (present(maskout)) then
961 0 : allocate(maskout(var_shape(1), var_shape(2)))
962 0 : maskout = abs(real(data, dp) - nodata_value) .gt. tiny(1.0_dp)
963 : end if
964 :
965 0 : call nc%close()
966 0 : end subroutine read_spatial_data_nc_i4
967 :
968 :
969 : ! NAME
970 : ! read_spatial_data_nc_dp
971 :
972 : ! PURPOSE
973 : !> \brief Read file from filepath. If there is a nc file with the correct name it is prefered over the asci file.
974 :
975 : !> \details TODO: add description
976 :
977 : ! HISTORY
978 : !> \authors Simon Lüdke
979 :
980 : !> \date June 2025
981 0 : subroutine read_spatial_data_nc_dp(ncname, data, maskout, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value, varName)
982 0 : use mo_netcdf, only : NcDataset, NcVariable
983 : use mo_message, only : error_message
984 : implicit none
985 :
986 : character(len=*), intent(in) :: ncname
987 : character(len=*), intent(in), optional :: varName
988 : real(dp), dimension(:, :), allocatable, intent(out) :: data
989 : logical, optional, allocatable, dimension(:, :), intent(out) :: maskout
990 : integer(i4), intent(out) :: ncols, nrows
991 : real(dp), intent(out) :: xllcorner, yllcorner, cellsize
992 : real(dp), intent(out) :: nodata_value
993 :
994 : type(NcDataset) :: nc
995 : type(NcVariable) :: var
996 0 : integer(i4), allocatable :: var_shape(:), start(:), cnt(:)
997 : character(len=256) :: varName_local
998 :
999 : ! Open NetCDF dataset
1000 0 : nc = NcDataset(ncname, "r")
1001 :
1002 0 : if (present(varName)) then
1003 0 : varName_local = trim(varName)
1004 : else
1005 0 : varName_local = determine_data_var_name(nc, ncname)
1006 : end if
1007 :
1008 : ! Extract header info
1009 0 : call get_header_info_from_nc(ncname, 0_i4, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value, var=varName_local)
1010 :
1011 : ! Retrieve variable
1012 0 : var = nc%getVariable(trim(varName_local))
1013 :
1014 : ! Determine shape and prepare read extents
1015 0 : var_shape = var%getShape()
1016 0 : allocate(start(size(var_shape)), source=1_i4)
1017 0 : allocate(cnt(size(var_shape)), source=1_i4)
1018 0 : cnt(1:2) = var_shape(1:2)
1019 :
1020 : ! Determine no-data value
1021 0 : if (var%hasAttribute("_FillValue")) then
1022 0 : call var%getAttribute("_FillValue", nodata_value)
1023 0 : else if (var%hasAttribute("missing_value")) then
1024 0 : call var%getAttribute("missing_value", nodata_value)
1025 : else
1026 0 : call error_message('***ERROR: read_nc_i4_data: missing _FillValue or missing_value attribute')
1027 : end if
1028 :
1029 : ! Allocate and read data
1030 0 : allocate(data(var_shape(1), var_shape(2)))
1031 0 : call var%getData(data, start=start, cnt=cnt)
1032 :
1033 0 : if (present(maskout)) then
1034 0 : allocate(maskout(var_shape(1), var_shape(2)))
1035 0 : maskout = abs(data - nodata_value) .gt. tiny(1.0_dp)
1036 : end if
1037 :
1038 0 : call nc%close()
1039 0 : end subroutine read_spatial_data_nc_dp
1040 :
1041 :
1042 :
1043 : ! NAME
1044 : ! read_spatial_data_nc_or_ascii_dp
1045 :
1046 : ! PURPOSE
1047 : !> \brief Read file from filepath. If there is a nc file with the correct name it is prefered over the asci file.
1048 :
1049 : !> \details TODO: add description
1050 :
1051 : ! HISTORY
1052 : !> \authors Simon Lüdke
1053 :
1054 : !> \date June 2025
1055 :
1056 66 : subroutine read_spatial_data_nc_or_ascii_dp(filepath, fileunit, header_ncols, header_nrows, header_xllcorner, &
1057 : header_yllcorner, header_cellsize, data, maskout, &
1058 : out_nCols, out_nRows, out_xllcorner, out_yllcorner, out_cellsize, out_nodata_value)
1059 0 : use mo_os, only : path_root, path_isfile
1060 : implicit none
1061 :
1062 : ! filename with location
1063 : character(len = *), intent(in) :: filepath
1064 :
1065 : ! unit for opening the file
1066 : integer(i4), intent(in) :: fileunit
1067 :
1068 : ! number of rows of data fields:
1069 : integer(i4), intent(in) :: header_nRows
1070 :
1071 : ! number of columns of data fields:
1072 : integer(i4), intent(in) :: header_nCols
1073 :
1074 : ! header read in lower left corner
1075 : real(dp), intent(in) :: header_xllcorner
1076 :
1077 : ! header read in lower left corner
1078 : real(dp), intent(in) :: header_yllcorner
1079 :
1080 : ! header read in cellsize
1081 : real(dp), intent(in) :: header_cellsize
1082 :
1083 : ! data
1084 : real(dp), dimension(:, :), allocatable, intent(out) :: data
1085 :
1086 : ! mask
1087 : logical, optional, dimension(:, :), allocatable, intent(out) :: maskout
1088 :
1089 : ! number of rows of data fields:
1090 : integer(i4), optional, intent(out) :: out_nRows
1091 :
1092 : ! number of columns of data fields:
1093 : integer(i4), optional, intent(out) :: out_nCols
1094 :
1095 : ! header read in lower left corner
1096 : real(dp), optional, intent(out) :: out_xllcorner
1097 :
1098 : ! header read in lower left corner
1099 : real(dp), optional, intent(out) :: out_yllcorner
1100 :
1101 : ! header read in cellsize
1102 : real(dp), optional, intent(out) :: out_cellsize
1103 :
1104 : !
1105 : real(dp), optional, intent(out) :: out_nodata_value
1106 :
1107 : ! file exists
1108 66 : real(dp) :: nodata_value
1109 : integer(i4) :: nrows, ncols
1110 66 : real(dp) :: xllcorner, yllcorner, cellsize
1111 66 : character(len=:), allocatable :: ncname
1112 :
1113 : ! Prefer explicit NetCDF path if it already ends with .nc, otherwise use sidecar .nc
1114 66 : if (len_trim(filepath) >= 3) then
1115 66 : if ((filepath(len_trim(filepath)-2:len_trim(filepath)) == '.nc') .or. (filepath(len_trim(filepath)-2:len_trim(filepath)) == '.NC')) then
1116 0 : ncname = trim(filepath)
1117 : else
1118 66 : ncname = path_root(filepath) // '.nc'
1119 : end if
1120 : else
1121 0 : ncname = path_root(filepath) // '.nc'
1122 : end if
1123 :
1124 : ! preferably use nc file if it exists alternatively the asci version
1125 : ! print *, "Check if ", ncname, " existis: ", path_isfile(ncname)
1126 132 : if (path_isfile(ncname)) then
1127 0 : print *, "Read nc file: ", ncname
1128 0 : call read_spatial_data_nc(ncname, data, maskout, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value)
1129 : else
1130 66 : print *, "Read ascii file: ", filepath
1131 66 : if ((header_ncols /= 0_i4) .and. (header_nrows /= 0_i4)) then
1132 66 : ncols = header_ncols
1133 66 : nrows = header_nrows
1134 66 : yllcorner = header_yllcorner
1135 66 : xllcorner = header_xllcorner
1136 66 : cellsize = header_cellsize
1137 : else
1138 : call read_header_ascii(filepath, fileunit, &
1139 : ncols, nrows, xllcorner, &
1140 0 : yllcorner, cellsize, nodata_value)
1141 :
1142 : end if
1143 : ! print *, "run ascii dp reading with ", filepath
1144 : call read_spatial_data_ascii(filepath, fileunit, ncols, nrows, xllcorner, &
1145 66 : yllcorner, cellsize, data, maskout)
1146 : end if
1147 : ! if header information is needed as output it is saved to the output variables
1148 66 : if (present(out_nCols)) &
1149 0 : out_nCols = ncols
1150 66 : if (present(out_nRows)) &
1151 0 : out_nRows = nrows
1152 66 : if (present(out_xllcorner)) &
1153 0 : out_xllcorner = xllcorner
1154 66 : if (present(out_yllcorner)) &
1155 0 : out_yllcorner = yllcorner
1156 66 : if (present(out_cellsize)) &
1157 0 : out_cellsize = cellsize
1158 66 : if (present(out_nodata_value)) &
1159 0 : out_nodata_value = nodata_value
1160 :
1161 66 : end subroutine read_spatial_data_nc_or_ascii_dp
1162 :
1163 : ! NAME
1164 : ! read_spatial_data_nc_or_ascii_i4
1165 :
1166 : ! PURPOSE
1167 : !> \brief Read file from filepath. If there is a nc file with the correct name it is prefered over the asci file.
1168 :
1169 : !> \details TODO: add description
1170 :
1171 : ! HISTORY
1172 : !> \authors Simon Lüdke
1173 :
1174 : !> \date June 2025
1175 :
1176 192 : subroutine read_spatial_data_nc_or_ascii_i4(filepath, fileunit, header_ncols, header_nrows, header_xllcorner, &
1177 : header_yllcorner, header_cellsize, data, maskout, &
1178 : out_nCols, out_nRows, out_xllcorner, out_yllcorner, out_cellsize, out_nodata_value)
1179 66 : use mo_os, only : path_root, path_isfile
1180 : implicit none
1181 :
1182 : ! filename with location
1183 : character(len = *), intent(in) :: filepath
1184 :
1185 : ! unit for opening the file
1186 : integer(i4), intent(in) :: fileunit
1187 :
1188 : ! number of rows of data fields:
1189 : integer(i4), intent(in) :: header_nRows
1190 :
1191 : ! number of columns of data fields:
1192 : integer(i4), intent(in) :: header_nCols
1193 :
1194 : ! header read in lower left corner
1195 : real(dp), intent(in) :: header_xllcorner
1196 :
1197 : ! header read in lower left corner
1198 : real(dp), intent(in) :: header_yllcorner
1199 :
1200 : ! header read in cellsize
1201 : real(dp), intent(in) :: header_cellsize
1202 :
1203 : ! data
1204 : integer(i4), dimension(:, :), allocatable, intent(out) :: data
1205 :
1206 : ! mask
1207 : logical, optional, dimension(:, :), allocatable, intent(out) :: maskout
1208 :
1209 : ! number of rows of data fields:
1210 : integer(i4), optional, intent(out) :: out_nRows
1211 :
1212 : ! number of columns of data fields:
1213 : integer(i4), optional, intent(out) :: out_nCols
1214 :
1215 : ! header read in lower left corner
1216 : real(dp), optional, intent(out) :: out_xllcorner
1217 :
1218 : ! header read in lower left corner
1219 : real(dp), optional, intent(out) :: out_yllcorner
1220 :
1221 : ! header read in cellsize
1222 : real(dp), optional, intent(out) :: out_cellsize
1223 :
1224 : real(dp), optional, intent(out) :: out_nodata_value
1225 :
1226 : ! file exists
1227 192 : real(dp) :: nodata_value
1228 : integer(i4) :: nrows, ncols
1229 192 : real(dp) :: xllcorner, yllcorner, cellsize
1230 192 : character(len=:), allocatable :: ncname
1231 :
1232 : ! Prefer explicit NetCDF path if it already ends with .nc, otherwise use sidecar .nc
1233 192 : if (len_trim(filepath) >= 3) then
1234 192 : if ((filepath(len_trim(filepath)-2:len_trim(filepath)) == '.nc') .or. (filepath(len_trim(filepath)-2:len_trim(filepath)) == '.NC')) then
1235 0 : ncname = trim(filepath)
1236 : else
1237 192 : ncname = path_root(filepath) // '.nc'
1238 : end if
1239 : else
1240 0 : ncname = path_root(filepath) // '.nc'
1241 : end if
1242 :
1243 : ! preferably use nc file if it exists alternatively the asci version
1244 : ! print *, "Check if ", ncname, " existis: ", path_isfile(ncname)
1245 384 : if (path_isfile(ncname)) then
1246 0 : print *, "Read nc file: ", ncname
1247 0 : call read_spatial_data_nc(ncname, data, maskout, ncols, nrows, xllcorner, yllcorner, cellsize, nodata_value)
1248 : else
1249 192 : print *, "Read ascii file: ", filepath
1250 192 : if ((header_ncols /= 0_i4) .and. (header_nrows /= 0_i4)) then
1251 192 : ncols = header_ncols
1252 192 : nrows = header_nrows
1253 192 : yllcorner = header_yllcorner
1254 192 : xllcorner = header_xllcorner
1255 192 : cellsize = header_cellsize
1256 : else
1257 : call read_header_ascii(filepath, fileunit, &
1258 : ncols, nrows, xllcorner, &
1259 0 : yllcorner, cellsize, nodata_value)
1260 :
1261 : end if
1262 : ! print *, "run ascii i4 reading with ", filepath, " as ", ncname, " does not exists"
1263 : call read_spatial_data_ascii(filepath, fileunit, ncols, nrows, xllcorner, &
1264 192 : yllcorner, cellsize, data, maskout)
1265 : end if
1266 : ! if header information is needed as output it is saved to the output variables
1267 192 : if (present(out_nCols)) &
1268 0 : out_nCols = ncols
1269 192 : if (present(out_nRows)) &
1270 0 : out_nRows = nrows
1271 192 : if (present(out_xllcorner)) &
1272 0 : out_xllcorner = xllcorner
1273 192 : if (present(out_yllcorner)) &
1274 0 : out_yllcorner = yllcorner
1275 192 : if (present(out_cellsize)) &
1276 0 : out_cellsize = cellsize
1277 192 : if (present(out_nodata_value)) &
1278 0 : out_nodata_value = nodata_value
1279 192 : end subroutine read_spatial_data_nc_or_ascii_i4
1280 :
1281 :
1282 0 : END MODULE mo_read_spatial_data
|
