libsim Versione 7.2.6
grid_id_class.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18#include "config.h"
19!> This module defines an abstract interface to different drivers for
20!! access to files containing gridded information. It defines a class
21!! associated to a file-like object \a grid_file_id, and a class
22!! associated to a grid-like object \a grid_id, extracted from a
23!! grid_file_id object. At the moment it has support for grib_api and
24!! gdal.
25!!
26!! Example of use (driver independent):
27!! \code
28!! ...
29!! USE grid_id_class
30!!
31!! TYPE(grid_file_id) :: ifile, ofile
32!! TYPE(grid_id) :: gaid
33!! CHARACTER(len=512) :: ifilename, ofilename
34!!
35!! ! open files
36!! ifile = grid_file_id_new(trim(ifilename),'r')
37!! ofile = grid_file_id_new(trim(ofilename),'w')
38!! IF (.NOT.c_e(ifile) .OR. .NOT.c_e(ofile)) STOP
39!!
40!! ! loop on all the messages in a file
41!! DO WHILE (.TRUE.)
42!! ! import from input file
43!! gaid = grid_id_new(ifile)
44!! IF (.NOT.c_e(gaid)) EXIT
45!!
46!! ! work on gaid here
47!! ...
48!! ! export to output file
49!! CALL export(gaid, ofile)
50!! CALL delete(gaid)
51!! ENDDO
52!!
53!! ! close the files
54!! CALL delete(ifile)
55!! CALL delete(ofile)
56!! ...
57!!
58!! \endcode
59!!
60!! \ingroup volgrid6d
61MODULE grid_id_class
62#ifdef HAVE_LIBGRIBAPI
63USE grib_api
64#endif
65#ifdef HAVE_LIBGDAL
66USE gdal
67#endif
74IMPLICIT NONE
75
76INTEGER,PARAMETER :: grid_id_no_driver = 0 !< constants to be used for associating an object to a driver: no type specified
77INTEGER,PARAMETER :: grid_id_grib_api = 1 !< type grib_api specified
78INTEGER,PARAMETER :: grid_id_gdal = 2 !< type gdal specified
79
80#if defined HAVE_LIBGRIBAPI
81INTEGER,PARAMETER :: grid_id_default = grid_id_grib_api !< default driver if none specified in constructor
82#elif defined HAVE_LIBGDAL
83INTEGER,PARAMETER :: grid_id_default = grid_id_gdal !< default driver if none specified in constructor
84#else
85INTEGER,PARAMETER :: grid_id_default = grid_id_no_driver !< default driver if none specified in constructor
86#endif
87
88CHARACTER(len=12),PARAMETER :: driverlist(0:2) = &
89 (/'no_driver ','grib_api ','gdal '/)
90
91#ifdef HAVE_LIBGDAL
92!> Derived type containing driver-specific options for gdal.
93!! It is initialised when opening the file/dataset with the
94!! constructor grid_file_id_new with the information provided in the
95!! driver string.
97 DOUBLE PRECISION :: xmin=dmiss !< bounding box of the area to import from the dataset
98 DOUBLE PRECISION :: ymin=dmiss !< bounding box of the area to import from the dataset
99 DOUBLE PRECISION :: xmax=dmiss !< bounding box of the area to import from the dataset
100 DOUBLE PRECISION :: ymax=dmiss !< bounding box of the area to import from the dataset
102#endif
103
104!> Derived type associated to a file-like object containing many
105!! blocks/messages/records/bands of gridded data.
106TYPE grid_file_id
107PRIVATE
108#ifdef HAVE_LIBGRIBAPI
109INTEGER :: gaid=imiss
110#endif
111#ifdef HAVE_LIBGDAL
112TYPE(gdaldataseth) :: gdalid
113INTEGER :: nlastband=0
114TYPE(gdal_file_id_options) :: gdal_options
115TYPE(grid_file_id),POINTER :: file_id_copy=>null()
116#endif
117INTEGER :: driver=grid_id_default
118END TYPE grid_file_id
119
120
121!> Derived type associated to a block/message/record/band of gridded
122!! data coming from a file-like object.
123TYPE grid_id
124PRIVATE
125INTEGER :: nodriverid=imiss
126#ifdef HAVE_LIBGRIBAPI
127INTEGER :: gaid=imiss
128#endif
129#ifdef HAVE_LIBGDAL
130TYPE(gdalrasterbandh) :: gdalid
131TYPE(grid_file_id),POINTER :: file_id=>null()
132#endif
133INTEGER :: driver=grid_id_default
134END TYPE grid_id
135
136!> Constructors for the corresponding classes in SUBROUTINE form.
137!! It is alternative to the *_new function constructors.
138INTERFACE init
139 MODULE PROCEDURE grid_file_id_init, grid_id_init
140END INTERFACE
141
142!> Destructors for the corresponding classes.
143INTERFACE delete
144 MODULE PROCEDURE grid_file_id_delete, grid_id_delete
145END INTERFACE
146
147!> Make a deep copy, if possible, of the grid identifier.
148INTERFACE copy
149 MODULE PROCEDURE grid_id_copy
150END INTERFACE
151
152!> Export a grid to a file
153INTERFACE export
154 MODULE PROCEDURE grid_id_export
155END INTERFACE
156
157!INTERFACE ASSIGNMENT(=)
158! MODULE PROCEDURE &
159!#ifdef HAVE_LIBGRIBAPI
160! grid_id_assign_int, &
161!#endif
162!#ifdef HAVE_LIBGDAL
163! grid_id_assign_dal, &
164!this%gdalid = that%gdalid
165!#endif
166! grid_id_assign_grid_id
167!END INTERFACE
168
169!> Check whether the corresponding object has been correctly associated.
170
171!> Set of LOGICAL functions to check whether a \a grid_file_id or a \a
172!! grid_id object have been correctly associated to a file or a grid.
173!! For a \a grid_file_id object it returns \a .FALSE. if the file has
174!! not been opened correctly or if the object has been initialized as
175!! empty. For a \a grid_id object it returns \a .FALSE. if the grid
176!! has not been correctly obtained from the file or if the object has
177!! been initialized as empty. They work both on scalars and 1-d array
178!! objects.
179!!
180!! \param this (TYPE(grid_file_id) or TYPE(grid_id)) object to be checked.
181INTERFACE c_e
182 MODULE PROCEDURE grid_id_c_e, grid_id_c_e_v, grid_file_id_c_e, grid_file_id_c_e_v
183END INTERFACE
184
185!> Display on standard output a description of the \a grid_id object
186!! provided.
187INTERFACE display
188 MODULE PROCEDURE grid_id_display
189END INTERFACE
190
191PRIVATE grid_file_id_delete, grid_id_delete, grid_id_copy, &
192 grid_id_c_e, grid_file_id_c_e, grid_id_c_e_v, grid_file_id_c_e_v, grid_id_display
193
194CONTAINS
195
196
197SUBROUTINE grid_file_id_init(this, filename, mode, driver, from_grid_id)
198TYPE(grid_file_id),INTENT(out) :: this ! object to initialise
199CHARACTER(len=*),INTENT(in) :: filename ! name of file containing gridded data, in the format [driver:]pathname
200CHARACTER(len=*),INTENT(in) :: mode ! access mode for file, 'r' or 'w'
201INTEGER,INTENT(in),OPTIONAL :: driver ! select the driver that will be associated to the grid_file_id created, use the constants \a grid_id_notype, \a grid_id_grib_api, \a grid_id_gdal
202TYPE(grid_id),INTENT(in),OPTIONAL :: from_grid_id ! select the driver as the one associated to the provided grid_id object
203
204this = grid_file_id_new(filename, mode, driver, from_grid_id)
205
206END SUBROUTINE grid_file_id_init
207
208
209!> Constructor for the \a grid_file_id class. It opens the associated
210!! file(s); the driver to be used for file access is selected
211!! according to the \a filename argument, to the optional argument \a
212!! driver, or to the optional argument \a from_grid_id, with
213!! increasing priority. If \a driver and \a from_grid_id are not
214!! provided and \a filename does not contain driver information, a
215!! default is chosen. If filename is an empty string or missing value,
216!! the object will be empty, the same will happen in case the file
217!! cannot be successfully opened. This condition can be tested with
218!! the function \a c_e() . The driver string provided with the
219!! filename can also contain driver-specific options separated by
220!! commas, e.g. \c 'gdal,8,44,10,46:globe.dat'.
221FUNCTION grid_file_id_new(filename, mode, driver, from_grid_id) RESULT(this)
222CHARACTER(len=*),INTENT(in) :: filename !< name of file containing gridded data, in the format [driver:]pathname
223CHARACTER(len=*),INTENT(in) :: mode !< access mode for file, 'r' or 'w'
224INTEGER,INTENT(in),OPTIONAL :: driver !< select the driver that will be associated to the grid_file_id created, use the constants \a grid_id_notype, \a grid_id_grib_api, \a grid_id_gdal
225TYPE(grid_id),INTENT(in),OPTIONAL :: from_grid_id !< select the driver as the one associated to the provided grid_id object
226TYPE(grid_file_id) :: this
227
228INTEGER :: n, ier, nf
229#ifdef HAVE_LIBGDAL
230INTEGER :: imode
231#endif
232TYPE(csv_record) :: driveropts
233CHARACTER(len=12) :: drivername
234
235#ifdef HAVE_LIBGDAL
236CALL gdalnullify(this%gdalid)
237#endif
238
239IF (filename == '' .OR. .NOT.c_e(filename)) RETURN
240
241n = index(filename,':')
242IF (n > 1) THEN ! override with driver from filename
243 CALL init(driveropts, filename(:n-1), nfield=nf)
244 CALL csv_record_getfield(driveropts, drivername)
245#ifdef HAVE_LIBGRIBAPI
246 IF (drivername == 'grib_api') THEN
247 this%driver = grid_id_grib_api
248 ENDIF
249#endif
250#ifdef HAVE_LIBGDAL
251 IF (drivername == 'gdal') THEN
252 IF (nf > 4) THEN
253 this%driver = grid_id_gdal
254 CALL csv_record_getfield(driveropts, this%gdal_options%xmin)
255 CALL csv_record_getfield(driveropts, this%gdal_options%ymin)
256 CALL csv_record_getfield(driveropts, this%gdal_options%xmax)
257 CALL csv_record_getfield(driveropts, this%gdal_options%ymax)
258 ! set extreme values if missing, magnitude defined empirically
259 ! because of integer overflow in fortrangis
260 IF (.NOT.c_e(this%gdal_options%xmin)) this%gdal_options%xmin = -1.0d6
261 IF (.NOT.c_e(this%gdal_options%ymin)) this%gdal_options%ymin = -1.0d6
262 IF (.NOT.c_e(this%gdal_options%xmax)) this%gdal_options%xmax = 1.0d6
263 IF (.NOT.c_e(this%gdal_options%ymax)) this%gdal_options%ymax = 1.0d6
264 ELSE
265 CALL l4f_log(l4f_error, 'gdal driver requires 4 extra arguments (bounding box)')
266 CALL raise_error()
267 ENDIF
268 ENDIF
269#endif
270 CALL delete(driveropts)
271ENDIF
272IF (PRESENT(driver)) THEN ! override with driver
273 this%driver = driver
274ENDIF
275IF (PRESENT(from_grid_id)) THEN ! override with from_grid_id
276 this%driver = from_grid_id%driver
277ENDIF
278
279#ifdef HAVE_LIBGRIBAPI
280IF (this%driver == grid_id_grib_api) THEN
281 CALL grib_open_file(this%gaid, filename(n+1:), trim(mode), ier)
282 IF (ier /= grib_success) this%gaid = imiss
283ENDIF
284#endif
285#ifdef HAVE_LIBGDAL
286IF (this%driver == grid_id_gdal) THEN
287 IF (mode(1:1) == 'w') THEN
288 imode = ga_update
289 ELSE
290 imode = ga_readonly
291 ENDIF
292 CALL gdalallregister()
293 this%gdalid = gdalopen(trim(filename(n+1:))//c_null_char, imode)
294! dirty trick, with gdal I have to keep a copy of the file_id, memory leak
295 ALLOCATE(this%file_id_copy)
296 this%file_id_copy = this
297ENDIF
298#endif
299
300END FUNCTION grid_file_id_new
301
302
303!> Count the number of block/message/record/band of gridded
304!! data in the file-like object provided. Returns 0 if the \a file_id
305!! object is empty or not corrctly associated to a file.
306FUNCTION grid_file_id_count(this) RESULT(count)
307TYPE(grid_file_id),INTENT(in) :: this !< file object to count
308INTEGER :: count
309
310INTEGER :: ier
312count = 0
313#ifdef HAVE_LIBGRIBAPI
314IF (this%driver == grid_id_grib_api) THEN
315 IF (c_e(this%gaid)) THEN
316 CALL grib_count_in_file(this%gaid, count, ier)
317 IF (ier /= grib_success) count = 0
318 ENDIF
319ENDIF
320#endif
321#ifdef HAVE_LIBGDAL
322IF (this%driver == grid_id_gdal) THEN
323 IF (gdalassociated(this%gdalid)) THEN
324 count = gdalgetrastercount(this%gdalid)
325 ENDIF
326ENDIF
327#endif
328
329END FUNCTION grid_file_id_count
330
332!> Destructor for the \a grid_file_id class. It closes the associated
333!! file(s) and releases the associated memory, but, in some drivers
334!! like grib_api, it may not release the memory associated to the \a
335!! grid_id objects read from that file, which continue their life in
336!! memory.
337SUBROUTINE grid_file_id_delete(this)
338TYPE(grid_file_id),INTENT(inout) :: this !< object to be deleted
339
340#ifdef HAVE_LIBGRIBAPI
341IF (this%driver == grid_id_grib_api) THEN
342 IF (c_e(this%gaid)) CALL grib_close_file(this%gaid)
343ENDIF
344this%gaid = imiss
345#endif
346#ifdef HAVE_LIBGDAL
347IF (this%driver == grid_id_gdal) THEN
348! dirty trick, with gdal I have to keep the file open
349! IF (gdalassociated(this%gdalid)) CALL gdalclose(this%gdalid)
350 this%nlastband = 0
351ENDIF
352CALL gdalnullify(this%gdalid)
353#endif
354
355this%driver = imiss
356
357END SUBROUTINE grid_file_id_delete
358
359
360! Function to check whether a \a grid_file_id object has been correctly associated
361! to the requested file. It returns \a .FALSE. if the file has not
362! been opened correctly or if the object has been initialized as empty.
363FUNCTION grid_file_id_c_e(this)
364TYPE(grid_file_id),INTENT(in) :: this ! object to be checked
365LOGICAL :: grid_file_id_c_e
366
367grid_file_id_c_e = .false.
368
369#ifdef HAVE_LIBGRIBAPI
370IF (this%driver == grid_id_grib_api) THEN
371 grid_file_id_c_e = c_e(this%gaid)
372ENDIF
373#endif
374#ifdef HAVE_LIBGDAL
375IF (this%driver == grid_id_gdal) THEN
376 grid_file_id_c_e = gdalassociated(this%gdalid)
377ENDIF
378#endif
379
380END FUNCTION grid_file_id_c_e
381
382
383! Function to check whether a \a grid_file_id object has been correctly associated
384! to the requested file. It returns \a .FALSE. if the file has not
385! been opened correctly or if the object has been initialized as empty.
386FUNCTION grid_file_id_c_e_v(this)
387TYPE(grid_file_id),INTENT(in) :: this(:) ! object to be checked
388LOGICAL :: grid_file_id_c_e_v(SIZE(this))
389
390INTEGER :: i
391
392DO i = 1, SIZE(this)
393 grid_file_id_c_e_v(i) = c_e(this(i))
394ENDDO
395
396END FUNCTION grid_file_id_c_e_v
397
398
399SUBROUTINE grid_id_init(this, from_grid_file_id, grib_api_template, grib_api_id)
400TYPE(grid_id),INTENT(out) :: this ! object to be initialized
401TYPE(grid_file_id),INTENT(inout),OPTIONAL :: from_grid_file_id ! file object from which grid object has to be created
402CHARACTER(len=*),INTENT(in),OPTIONAL :: grib_api_template ! grib_api template file from which grid_object has to be created
403INTEGER,INTENT(in),OPTIONAL :: grib_api_id ! grib_api id obtained directly from a \a grib_get subroutine call
404
405this = grid_id_new(from_grid_file_id, grib_api_template, grib_api_id)
406
407END SUBROUTINE grid_id_init
408
410!> Constructor for the \a grid_id class. It gets the next grid (grib
411!! message or raster band) from the file_id provided. If the file
412!! associated to the file_id provided contains no more grids, or if
413!! the argument \a file_id is not provided, an empty object is
414!! created; this condition can be tested with the function c_e().
415!! Alternative ways to define the object (to be used in rare cases)
416!! are through a grib_api template file name (\a grib_api_template
417!! argument) or through a grib_api integer id obtained directly from
418!! grib_api calls (\a grib_api_id argument).
419FUNCTION grid_id_new(from_grid_file_id, grib_api_template, grib_api_id, &
420 no_driver_id) RESULT(this)
421TYPE(grid_file_id),INTENT(inout),OPTIONAL,TARGET :: from_grid_file_id !< file object from which grid object has to be created
422CHARACTER(len=*),INTENT(in),OPTIONAL :: grib_api_template !< grib_api template file from which grid_object has to be created
423INTEGER,INTENT(in),OPTIONAL :: grib_api_id !< grib_api id obtained directly from a \a grib_get subroutine call
424INTEGER,INTENT(in),OPTIONAL :: no_driver_id
425TYPE(grid_id) :: this
426
427INTEGER :: ier
428
429#ifdef HAVE_LIBGDAL
430CALL gdalnullify(this%gdalid)
431#endif
432
433IF (PRESENT(from_grid_file_id)) THEN
434 this%driver = from_grid_file_id%driver ! take driver from file_id
435
436#ifdef HAVE_LIBGRIBAPI
437 IF (this%driver == grid_id_grib_api) THEN
438 IF (c_e(from_grid_file_id%gaid)) THEN
439 CALL grib_new_from_file(from_grid_file_id%gaid, this%gaid, ier)
440 IF (ier /= grib_success) this%gaid = imiss
441 ENDIF
442 ENDIF
443#endif
444#ifdef HAVE_LIBGDAL
445 IF (this%driver == grid_id_gdal) THEN
446 IF (gdalassociated(from_grid_file_id%gdalid) .AND. &
447 ASSOCIATED(from_grid_file_id%file_id_copy)) THEN
448 IF (from_grid_file_id%nlastband < &
449 gdalgetrastercount(from_grid_file_id%gdalid)) THEN ! anything to read?
450 from_grid_file_id%nlastband = from_grid_file_id%nlastband + 1
451 this%gdalid = &
452 gdalgetrasterband(from_grid_file_id%gdalid, from_grid_file_id%nlastband)
453 this%file_id => from_grid_file_id%file_id_copy ! for gdal remember copy of file_id
454
455 ENDIF
456 ENDIF
457 ENDIF
458#endif
459
460#ifdef HAVE_LIBGRIBAPI
461ELSE IF (PRESENT(grib_api_template)) THEN
462 this%driver = grid_id_grib_api
463 CALL grib_new_from_samples(this%gaid, grib_api_template, ier)
464 IF (ier /= grib_success) this%gaid = imiss
465ELSE IF (PRESENT(grib_api_id)) THEN
466 this%driver = grid_id_grib_api
467 this%gaid = grib_api_id
468#endif
469ELSE IF (PRESENT(no_driver_id)) THEN
470 this%driver = grid_id_no_driver
471 this%nodriverid = no_driver_id
472ENDIF
473
474END FUNCTION grid_id_new
475
476
477!> Destructor for the \a grid_id class. It releases the memory associated with
478!! the grid descriptor identifier. In grib_api this is necessary and
479!! can be made also after closing the corresponding \a grid_file_id
480!! object; while for gdal this is a no-operation.
481SUBROUTINE grid_id_delete(this)
482TYPE(grid_id),INTENT(inout) :: this !< object to be deleted
483
484this%nodriverid = imiss
485#ifdef HAVE_LIBGRIBAPI
486IF (this%driver == grid_id_grib_api) THEN
487 IF (c_e(this%gaid)) CALL grib_release(this%gaid)
488ENDIF
489this%gaid = imiss
490#endif
491#ifdef HAVE_LIBGDAL
492CALL gdalnullify(this%gdalid)
493NULLIFY(this%file_id)
494#endif
495
496this%driver = imiss
497
498END SUBROUTINE grid_id_delete
499
500
501!> Check whether the grid_id object is readonly (.TRUE.) or allows
502!! writing bands (.FALSE.)
503FUNCTION grid_id_readonly(this) RESULT(readonly)
504TYPE(grid_id),INTENT(in) :: this !< object to test
505LOGICAL :: readonly
506
507readonly = this%driver /= grid_id_grib_api
508
509END FUNCTION grid_id_readonly
510
511
512!> Performs a "deep" copy of the \a grid_id object when possible.
513!! For grib_api this clones the grid_id generating a new independent
514!! object in memory, which can be manipulated without affecting the
515!! original one. The \a grid_id object \a that does not need to be
516!! initialized before the call.
517SUBROUTINE grid_id_copy(this, that)
518TYPE(grid_id),INTENT(in) :: this !< source object
519TYPE(grid_id),INTENT(out) :: that !< destination object, it must not be initialized
520
521that = this ! start with a shallow copy
522
523#ifdef HAVE_LIBGRIBAPI
524IF (this%driver == grid_id_grib_api) THEN
525 IF (c_e(this%gaid)) THEN
526 that%gaid = -1
527 CALL grib_clone(this%gaid, that%gaid)
528 ENDIF
529ENDIF
530#endif
531#ifdef HAVE_LIBGDAL
532IF (this%driver == grid_id_gdal) THEN
533 IF (c_e(this)) THEN
534! that = grid_id_new(no_driver_id=1)
535! that%gdalid = this%gdalid ! better idea?
536! that%file_id => this%file_id
537 ENDIF
538ENDIF
539#endif
540
541END SUBROUTINE grid_id_copy
542
543
544!> Export a \a grid_id object \a this to the file indicated by a
545!! \a grid_file_id object. Both grid_id and grid_file_id objects must
546!! be related to the same driver (e.g. grib_api or gdal).
547SUBROUTINE grid_id_export(this, file_id)
548TYPE(grid_id),INTENT(inout) :: this
549TYPE(grid_file_id),INTENT(in) :: file_id !< file object to which grid has to be exported
550
551INTEGER :: ier
552
553IF (c_e(this) .AND. c_e(file_id)) THEN
554#ifdef HAVE_LIBGRIBAPI
555 IF (this%driver == grid_id_grib_api .AND. file_id%driver == grid_id_grib_api) &
556 CALL grib_write(this%gaid, file_id%gaid, ier) ! log ier?
557#endif
558ENDIF
559#ifdef HAVE_LIBGDAL
560IF (this%driver == grid_id_gdal .AND. file_id%driver == grid_id_gdal) THEN
561 ! not implemented, log?
562ENDIF
563#endif
564
565END SUBROUTINE grid_id_export
566
567
568! Function to check whether a \a _file_id object has been correctly associated
569! to a grid. It returns \a .FALSE. if the grid has not been correctly
570! obtained from the file or if the object has been initialized as
571! empty.
572FUNCTION grid_id_c_e(this)
573TYPE(grid_id),INTENT(in) :: this ! object to be checked
574LOGICAL :: grid_id_c_e
575
576grid_id_c_e = .false.
577
578#ifdef HAVE_LIBGRIBAPI
579IF (this%driver == grid_id_grib_api) THEN
580 grid_id_c_e = c_e(this%gaid)
581ENDIF
582#endif
583#ifdef HAVE_LIBGDAL
584IF (this%driver == grid_id_gdal) THEN
585 grid_id_c_e = gdalassociated(this%gdalid)
586ENDIF
587#endif
588IF (this%driver == grid_id_no_driver) THEN
589 grid_id_c_e = c_e(this%nodriverid)
590ENDIF
591
592END FUNCTION grid_id_c_e
593
594
595! Function to check whether a \a _file_id object has been correctly associated
596! to a grid. It returns \a .FALSE. if the grid has not been correctly
597! obtained from the file or if the object has been initialized as
598! empty.
599FUNCTION grid_id_c_e_v(this)
600TYPE(grid_id),INTENT(in) :: this(:) ! object to be checked
601LOGICAL :: grid_id_c_e_v(SIZE(this))
602
603INTEGER :: i
604
605DO i = 1, SIZE(this)
606 grid_id_c_e_v(i) = c_e(this(i))
607ENDDO
608
609END FUNCTION grid_id_c_e_v
610
611
612!> Return a character representation of the driver associated
613!! with the object \a this. The result is a string with constant
614!! length and blank-padded at the end, representing the name of the
615!! driver; 'none' is retured if the object is empty.
616FUNCTION grid_file_id_get_driver(this) RESULT(driver)
617TYPE(grid_file_id),INTENT(in) :: this
618CHARACTER(len=LEN(driverlist)) :: driver
619
620IF (this%driver > 0 .AND. this%driver <= SIZE(driverlist)) THEN
621 driver = driverlist(this%driver)
622ELSE
623 driver = driverlist(0)
624ENDIF
625
626END FUNCTION grid_file_id_get_driver
627
628
629!> Return a character representation of the driver associated
630!! with the object \a this. The result is a string with constant
631!! length and blank-padded at the end, representing the name of the
632!! driver; 'none' is retured if the object is empty.
633FUNCTION grid_id_get_driver(this) RESULT(driver)
634TYPE(grid_id),INTENT(in) :: this
635CHARACTER(len=LEN(driverlist)) :: driver
636
637IF (this%driver > 0 .AND. this%driver <= SIZE(driverlist)) THEN
638 driver = driverlist(this%driver)
639ELSE
640 driver = driverlist(0)
641ENDIF
642
643END FUNCTION grid_id_get_driver
644
645
646!> Display on standard output a description of the \a grid_id object
647!! provided. Also the grib key names and values are printed; the set
648!! of keys returned can be controlled with the input variable
649!! namespace. Available namespaces are \c ls, to get the same default
650!! keys as the \a grib_ls command, and \c mars to get the keys used by
651!! mars.
652SUBROUTINE grid_id_display(this, namespace)
653TYPE(grid_id),INTENT(in) :: this !< object to be displayed
654CHARACTER(len=*),OPTIONAL :: namespace !< grib_api namespace of the keys to search for, all the keys if empty, default \c ls
655
656INTEGER :: kiter, iret
657CHARACTER(len=255) :: key, value, lnamespace
658
659
660#ifdef HAVE_LIBGRIBAPI
661IF (this%driver == grid_id_grib_api) THEN
662
663 lnamespace = optio_c(namespace,255)
664 IF (.NOT.c_e(lnamespace))THEN
665 lnamespace = "ls"
666 ENDIF
667
668 print*,"GRIB_API namespace:",trim(lnamespace)
670 CALL grib_keys_iterator_new(this%gaid, kiter, namespace=trim(lnamespace))
671
672 iter: DO
673 CALL grib_keys_iterator_next(kiter, iret)
674
675 IF (iret /= 1) THEN
676 EXIT iter
677 END IF
678
679 CALL grib_keys_iterator_get_name(kiter, key)
680!<\todo bug in grib_api, segmentation fault with computeStatistics key
681 IF (key == 'computeStatistics') cycle
682
683 CALL grib_get(this%gaid, trim(key), value, iret)
684 IF (iret == 0)THEN
685 print*, trim(key)//' = '//trim(VALUE)
686 ELSE
687 print*, trim(key)//' = '//"KEY NOT FOUND, namespace :"//trim(lnamespace)//" ( bug ? )"
688 ENDIF
689 ENDDO iter
690
691 CALL grib_keys_iterator_delete(kiter)
692
693ENDIF
694
695#endif
696
697END SUBROUTINE grid_id_display
698
699
700#ifdef HAVE_LIBGRIBAPI
701!> Returns the original grib_api id associated with the object
702!! provided, i.e. the unit associated by grib_api to the file.
703FUNCTION grid_file_id_get_gaid(this) RESULT(gaid)
704TYPE(grid_file_id),INTENT(in) :: this !< object to query
705INTEGER :: gaid
706gaid = this%gaid
707END FUNCTION grid_file_id_get_gaid
708
709!> Returns the original grib_api id associated with the object
710!! provided, i.e. the grib_api id associated to the grid.
711FUNCTION grid_id_get_gaid(this) RESULT(gaid)
712TYPE(grid_id),INTENT(in) :: this !< object to query
713INTEGER :: gaid
714gaid = this%gaid
715END FUNCTION grid_id_get_gaid
716#endif
717
718
719#ifdef HAVE_LIBGDAL
720!> Returns the original gdal Fortran object associated with the object
721!! provided, i.e. the dataset pointer.
722FUNCTION grid_file_id_get_gdalid(this) RESULT(gdalid)
723TYPE(grid_file_id),INTENT(in) :: this !< object to query
724TYPE(gdaldataseth) :: gdalid
725gdalid = this%gdalid
726END FUNCTION grid_file_id_get_gdalid
727
728!> Returns the original gdal Fortran object associated with the object
729!! provided, i.e. the raster band pointer.
730FUNCTION grid_id_get_gdalid(this) RESULT(gdalid)
731TYPE(grid_id),INTENT(in) :: this !< object to query
732TYPE(gdalrasterbandh) :: gdalid
733gdalid = this%gdalid
734END FUNCTION grid_id_get_gdalid
736!> Returns an object with driver-specific options associated to the grid_id
737!! provided, available only for gdal.
738FUNCTION grid_id_get_gdal_options(this) RESULT(gdal_options)
739TYPE(grid_id),INTENT(in) :: this !< object to query
740TYPE(gdal_file_id_options) :: gdal_options
741
742TYPE(gdal_file_id_options) :: gdal_options_local
743
744IF (ASSOCIATED(this%file_id)) THEN
745 gdal_options = this%file_id%gdal_options
746ELSE
747 gdal_options = gdal_options_local ! empty object
748ENDIF
749
750END FUNCTION grid_id_get_gdal_options
751#endif
752
753
754!> Decode and return the data array from a grid_id object.
755!! The output array \a field must have a size matching the size of the
756!! encoded data.
757SUBROUTINE grid_id_decode_data(this, field)
758TYPE(grid_id),INTENT(in) :: this
759REAL,INTENT(out) :: field(:,:)
760
761LOGICAL :: done
762
763done = .false.
764#ifdef HAVE_LIBGRIBAPI
765IF (c_e(this%gaid)) THEN
766 CALL grid_id_decode_data_gribapi(this%gaid, field)
767 done = .true.
768ENDIF
769#endif
770#ifdef HAVE_LIBGDAL
771! subarea?
772IF (gdalassociated(this%gdalid)) THEN
773 CALL grid_id_decode_data_gdal(this%gdalid, field, this%file_id%gdal_options)
774 done = .true.
775ENDIF
776#endif
777IF (.NOT.done) field(:,:) = rmiss
778
779END SUBROUTINE grid_id_decode_data
780
781
782!> Encode a data array into a grid_id object.
783!! The input array \a field must have a size matching the size of the
784!! dataset.
785SUBROUTINE grid_id_encode_data(this, field)
786TYPE(grid_id),INTENT(inout) :: this !< gridinfo object
787REAL,intent(in) :: field(:,:) !< data array to be encoded
788
789#ifdef HAVE_LIBGRIBAPI
790IF (this%driver == grid_id_grib_api) THEN
791
792! call display(this,"parameter")
793! print *,minval(field,c_e(field)),maxval(field,c_e(field))
794! print *,"-----------------------"
795
796 IF (c_e(this%gaid)) CALL grid_id_encode_data_gribapi(this%gaid, field)
797ENDIF
798#endif
799#ifdef HAVE_LIBGDAL
800IF (this%driver == grid_id_gdal) THEN
801!gdalid = grid_id_get_gdalid(this%gaid)
802 CALL l4f_log(l4f_warn,"export to gdal not implemented" )
803! subarea?
804ENDIF
805#endif
806
807END SUBROUTINE grid_id_encode_data
808
809
810#ifdef HAVE_LIBGRIBAPI
811SUBROUTINE grid_id_decode_data_gribapi(gaid, field)
812INTEGER,INTENT(in) :: gaid ! grib_api id
813REAL,INTENT(out) :: field(:,:) ! array of decoded values
814
815INTEGER :: editionnumber
816INTEGER :: alternativerowscanning, &
817 iscansnegatively, jscanspositively, jpointsareconsecutive
818INTEGER :: numberofvalues,numberofpoints
819REAL :: vector(size(field))
820INTEGER :: x1, x2, xs, y1, y2, ys, ord(2), ierr
822
823call grib_get(gaid,'GRIBEditionNumber',editionnumber)
824
825if (editionnumber == 2) then
826
827 CALL grib_get(gaid,'alternativeRowScanning',alternativerowscanning,ierr)
828 IF (ierr == grib_success .AND. alternativerowscanning /= 0) THEN
829 CALL l4f_log(l4f_error, "grib_api alternativeRowScanning not supported: " &
830 //t2c(alternativerowscanning))
831 CALL raise_error()
832 field(:,:) = rmiss
833 RETURN
834 ENDIF
835
836else if (editionnumber /= 1) then
837
838 CALL l4f_log(l4f_error, &
839 "grib_api GribEditionNumber not supported: "//t2c(editionnumber))
840 call raise_error()
841 field(:,:) = rmiss
842 RETURN
843
844end if
845
846CALL grib_get(gaid,'iScansNegatively',iscansnegatively,ierr)
847IF (ierr /= grib_success) iscansnegatively=0
848CALL grib_get(gaid,'jScansPositively',jscanspositively,ierr)
849IF (ierr /= grib_success) jscanspositively=1
850CALL grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive,ierr)
851IF (ierr /= grib_success) jpointsareconsecutive=0
852
853call grib_get(gaid,'numberOfPoints',numberofpoints)
854call grib_get(gaid,'numberOfValues',numberofvalues)
855
856IF (numberofpoints /= SIZE(field)) THEN
857 CALL l4f_log(l4f_error, 'grid_id_decode_data_gribapi numberOfPoints and grid size different')
858 CALL l4f_log(l4f_error, 'grid_id_decode_data_gribapi numberOfPoints: ' &
859 //t2c(numberofpoints)//', nx,ny: '&
860 //t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
861 CALL raise_error()
862 field(:,:) = rmiss
863 RETURN
864ENDIF
865
866! get data values
867#ifdef DEBUG
868call l4f_log(l4f_info,'grib_api number of values: '//to_char(numberofvalues))
869call l4f_log(l4f_info,'grib_api number of points: '//to_char(numberofpoints))
870#endif
871
872CALL grib_set(gaid,'missingValue',rmiss)
873CALL grib_get(gaid,'values',vector)
874! suspect bug in grib_api, when all field is missing it is set to zero
875IF (numberofvalues == 0) vector = rmiss
876
877#ifdef DEBUG
878CALL l4f_log(l4f_debug, 'grib_api, decoded field in interval: '// &
879 t2c(minval(vector,mask=c_e(vector)))//' '//t2c(maxval(vector,mask=c_e(vector))))
880CALL l4f_log(l4f_debug, 'grib_api, decoded field with number of missing: '// &
881 t2c(count(.NOT.c_e(vector))))
882#endif
883
884IF (numberofvalues /= count(c_e(vector))) THEN
885 CALL l4f_log(l4f_warn, 'grid_id_decode_data_gribapi numberOfValues and valid data count different')
886 CALL l4f_log(l4f_warn, 'grid_id_decode_data_gribapi numberOfValues: ' &
887 //t2c(numberofvalues)//', valid data: '//t2c(count(c_e(vector))))
888! CALL raise_warning()
889ENDIF
890
891! Transfer data field changing scanning mode to 64
892IF (iscansnegatively == 0) THEN
893 x1 = 1
894 x2 = SIZE(field,1)
895 xs = 1
896ELSE
897 x1 = SIZE(field,1)
898 x2 = 1
899 xs = -1
900ENDIF
901IF (jscanspositively == 0) THEN
902 y1 = SIZE(field,2)
903 y2 = 1
904 ys = -1
905ELSE
906 y1 = 1
907 y2 = SIZE(field,2)
908 ys = 1
909ENDIF
911IF ( jpointsareconsecutive == 0) THEN
912 ord = (/1,2/)
913ELSE
914 ord = (/2,1/)
915ENDIF
916
917field(x1:x2:xs,y1:y2:ys) = reshape(vector, &
918 (/SIZE(field,1),SIZE(field,2)/), order=ord)
919
920END SUBROUTINE grid_id_decode_data_gribapi
921
922
923SUBROUTINE grid_id_encode_data_gribapi(gaid, field)
924INTEGER,INTENT(in) :: gaid ! grib_api id
925REAL,intent(in) :: field(:,:) ! data array to be encoded
927INTEGER :: editionnumber
928INTEGER :: alternativerowscanning, iscansnegatively, &
929 jscanspositively, jpointsareconsecutive
930INTEGER :: x1, x2, xs, y1, y2, ys, ierr
931
932call grib_get(gaid,'GRIBEditionNumber',editionnumber)
933
934if (editionnumber == 2) then
935
936 CALL grib_get(gaid,'alternativeRowScanning',alternativerowscanning,ierr)
937 IF (ierr == grib_success .AND. alternativerowscanning /= 0)THEN
938 CALL l4f_log(l4f_error, "grib_api alternativeRowScanning not supported: " &
939 //trim(to_char(alternativerowscanning)))
940 CALL raise_error()
941 RETURN
942 ENDIF
943
944else if( editionnumber /= 1) then
946 call l4f_log(l4f_error, &
947 "grib_api GribEditionNumber not supported: "//t2c(editionnumber))
948 call raise_error()
949 RETURN
950
951end if
952
953CALL grib_get(gaid,'iScansNegatively',iscansnegatively,ierr)
954IF (ierr /= grib_success) iscansnegatively=0
955CALL grib_get(gaid,'jScansPositively',jscanspositively,ierr)
956IF (ierr /= grib_success) jscanspositively=1
957CALL grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive,ierr)
958IF (ierr /= grib_success) jpointsareconsecutive=0
959
960! these grib_sets are alredy done in gridinfo_export, but it seems
961! that it is necessary to repeat them here, they can fail with
962! unstructured grids, thus ierr
963#ifdef DEBUG
964CALL l4f_log(l4f_debug, 'grib_api, Ni,Nj:'//t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
965#endif
966CALL grib_set(gaid,'Ni',SIZE(field,1), ierr)
967CALL grib_set(gaid,'Nj',SIZE(field,2), ierr)
968
969! Transfer data field changing scanning mode from 64
970IF (iscansnegatively == 0) THEN
971 x1 = 1
972 x2 = SIZE(field,1)
973 xs = 1
974ELSE
975 x1 = SIZE(field,1)
976 x2 = 1
977 xs = -1
978ENDIF
979IF (jscanspositively == 0) THEN
980 y1 = SIZE(field,2)
981 y2 = 1
982 ys = -1
983ELSE
984 y1 = 1
985 y2 = SIZE(field,2)
986 ys = 1
987ENDIF
988
989
990IF (any(field == rmiss)) THEN
991
992 CALL grib_set(gaid,'missingValue',rmiss)
993 IF (editionnumber == 1) THEN
994! enable bitmap in grib1
995! grib_api 1.9.9 went into an infinite loop with second order packing here
996! CALL grib_set(gaid,'packingType','grid_simple')
997! now it's fixed, leaving second order if present
998 CALL grib_set(gaid,"bitmapPresent",1)
999 ELSE
1000! enable bitmap in grib2
1001 CALL grib_set(gaid,"bitMapIndicator",0)
1002 ENDIF
1003
1004ELSE
1005
1006 IF (editionnumber == 1) THEN
1007! disable bitmap in grib1
1008 CALL grib_set(gaid,"bitmapPresent",0)
1009 ELSE
1010! disable bitmap in grib2
1011 CALL grib_set(gaid,"bitMapIndicator",255)
1012 ENDIF
1013
1014ENDIF
1015
1016#ifdef DEBUG
1017CALL l4f_log(l4f_debug, 'grib_api, coding field in interval: '// &
1018 t2c(minval(field,mask=c_e(field)))//' '//t2c(maxval(field,mask=c_e(field))))
1019CALL l4f_log(l4f_debug, 'grib_api, coding field with number of missing: '// &
1020 t2c(count(.NOT.c_e(field))))
1021CALL l4f_log(l4f_debug, 'grib_api, sizex:'//t2c(x1)//','//t2c(x2)//','//t2c(xs))
1022CALL l4f_log(l4f_debug, 'grib_api, sizey:'//t2c(y1)//','//t2c(y2)//','//t2c(ys))
1023#endif
1024IF (jpointsareconsecutive == 0) THEN
1025 CALL grib_set(gaid,'values', reshape(field(x1:x2:xs,y1:y2:ys), &
1026 (/SIZE(field)/)))
1027ELSE
1028 CALL grib_set(gaid,'values', reshape(transpose(field(x1:x2:xs,y1:y2:ys)), &
1029 (/SIZE(field)/)))
1030ENDIF
1031
1032END SUBROUTINE grid_id_encode_data_gribapi
1033#endif
1034
1035
1036#ifdef HAVE_LIBGDAL
1037SUBROUTINE grid_id_decode_data_gdal(gdalid, field, gdal_options)
1038#ifdef F2003_FULL_FEATURES
1039USE ieee_arithmetic
1040#endif
1041TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal id
1042REAL,INTENT(out) :: field(:,:) ! array of decoded values
1043TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
1044
1045TYPE(gdaldataseth) :: hds
1046REAL(kind=c_double) :: geotrans(6), dummy1, dummy2, dummy3, dummy4
1047REAL :: gdalmiss
1048REAL,ALLOCATABLE :: buffer(:,:)
1049INTEGER :: ix1, iy1, ix2, iy2, ixs, iys, ord(2), ier
1050INTEGER(kind=c_int) :: nrx, nry
1051LOGICAL :: must_trans
1052
1053
1054hds = gdalgetbanddataset(gdalid) ! go back to dataset
1055ier = gdalgetgeotransform(hds, geotrans)
1056
1057IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN
1058! transformation is diagonal, no transposing
1059 IF (geotrans(2) > 0.0_c_double) THEN
1060 ix1 = 1
1061 ix2 = SIZE(field,1)
1062 ixs = 1
1063 ELSE
1064 ix1 = SIZE(field,1)
1065 ix2 = 1
1066 ixs = -1
1067 ENDIF
1068 IF (geotrans(6) > 0.0_c_double) THEN
1069 iy1 = 1
1070 iy2 = SIZE(field,2)
1071 iys = 1
1072 ELSE
1073 iy1 = SIZE(field,2)
1074 iy2 = 1
1075 iys = -1
1076 ENDIF
1077 nrx = SIZE(field,1)
1078 nry = SIZE(field,2)
1079 ord = (/1,2/)
1080 must_trans = .false.
1081! ALLOCATE(buffer(nrx,nry))
1082
1083ELSE IF (geotrans(2) == 0.0_c_double .AND. geotrans(6) == 0.0_c_double) THEN
1084! transformation is anti-diagonal, transposing required this should not happen
1085 IF (geotrans(3) > 0.0_c_double) THEN
1086 ix1 = 1
1087 ix2 = SIZE(field,1)
1088 ixs = 1
1089 ELSE
1090 ix1 = SIZE(field,1)
1091 ix2 = 1
1092 ixs = -1
1093 ENDIF
1094 IF (geotrans(5) > 0.0_c_double) THEN
1095 iy1 = 1
1096 iy2 = SIZE(field,2)
1097 iys = 1
1098 ELSE
1099 iy1 = SIZE(field,2)
1100 iy2 = 1
1101 iys = -1
1102 ENDIF
1103 nrx = SIZE(field,2)
1104 nry = SIZE(field,1)
1105 ord = (/2,1/)
1106 must_trans = .true.
1107! ALLOCATE(buffer(nry,nrx))
1108
1109ELSE ! transformation is a rotation, not supported
1110 CALL l4f_log(l4f_error, 'gdal geotransform is a generic rotation, not supported')
1111 CALL raise_error()
1112 field(:,:) = rmiss
1113 RETURN
1114ENDIF
1115
1116! read data from file
1117CALL gdalrastersimpleread_f(gdalid, gdal_options%xmin, gdal_options%ymin, &
1118 gdal_options%xmax, gdal_options%ymax, buffer, dummy1, dummy2, dummy3, dummy4)
1119
1120IF (.NOT.ALLOCATED(buffer)) THEN ! error in read
1121 CALL l4f_log(l4f_error, 'gdal error in reading with gdal driver')
1122 CALL raise_error()
1123 field(:,:) = rmiss
1124 RETURN
1125ENDIF
1126
1127IF (SIZE(buffer) /= (SIZE(field)))THEN
1128 CALL l4f_log(l4f_error, 'gdal raster band and gridinfo size different')
1129 CALL l4f_log(l4f_error, 'gdal rasterband: ' &
1130 //t2c(SIZE(buffer,1))//'X'//t2c(SIZE(buffer,2))//', nx,ny:' &
1131 //t2c(SIZE(field,ord(1)))//'X'//t2c(SIZE(field,ord(2))))
1132 CALL raise_error()
1133 field(:,:) = rmiss
1134 RETURN
1135ENDIF
1136
1137#ifdef F2003_FULL_FEATURES
1138! aparently gdal datasets may contain NaN
1139WHERE(ieee_is_nan(buffer))
1140 buffer = rmiss
1141END WHERE
1142#else
1143WHERE(buffer /= buffer)
1144 buffer = rmiss
1145END WHERE
1146#endif
1147
1148! set missing value if necessary
1149gdalmiss = real(gdalgetrasternodatavalue(gdalid, ier))
1150IF (ier /= 0) THEN ! success -> there are missing values
1151#ifdef DEBUG
1152 CALL l4f_log(l4f_info, 'gdal missing data value: '//trim(to_char(gdalmiss)))
1153#endif
1154 WHERE(buffer(:,:) == gdalmiss)
1155 buffer(:,:) = rmiss
1156 END WHERE
1157ELSE
1158#ifdef DEBUG
1159 CALL l4f_log(l4f_info, 'gdal no missing data found in band')
1160#endif
1161ENDIF
1162
1163! reshape the field
1164IF (must_trans) THEN
1165 field(ix1:ix2:ixs,iy1:iy2:iys) = transpose(buffer)
1166ELSE
1167 field(ix1:ix2:ixs,iy1:iy2:iys) = buffer(:,:)
1168ENDIF
1169
1170
1171END SUBROUTINE grid_id_decode_data_gdal
1172#endif
1173
1174END MODULE grid_id_class
Set of functions that return a trimmed CHARACTER representation of the input variable.
Set of functions that return a CHARACTER representation of the input variable.
Methods for successively obtaining the fields of a csv_record object.
Check whether the corresponding object has been correctly associated.
Make a deep copy, if possible, of the grid identifier.
Destructors for the corresponding classes.
Display on standard output a description of the grid_id object provided.
Export a grid to a file.
Constructors for the corresponding classes in SUBROUTINE form.
Index method.
Utilities for CHARACTER variables.
Gestione degli errori.
Utilities for managing files.
This module defines an abstract interface to different drivers for access to files containing gridded...
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Class for interpreting the records of a csv file.
Derived type containing driver-specific options for gdal.
Derived type associated to a file-like object containing many blocks/messages/records/bands of gridde...
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...

Generated with Doxygen.