libsim  Versione 7.1.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 
61 MODULE grid_id_class
62 #ifdef HAVE_LIBGRIBAPI
63 USE grib_api
64 #endif
65 #ifdef HAVE_LIBGDAL
66 USE gdal
67 #endif
72 USE log4fortran
73 USE err_handling
74 IMPLICIT NONE
75 
76 INTEGER,PARAMETER :: grid_id_no_driver = 0
77 INTEGER,PARAMETER :: grid_id_grib_api = 1
78 INTEGER,PARAMETER :: grid_id_gdal = 2
79 
80 #if defined HAVE_LIBGRIBAPI
81 INTEGER,PARAMETER :: grid_id_default = grid_id_grib_api
82 #elif defined HAVE_LIBGDAL
83 INTEGER,PARAMETER :: grid_id_default = grid_id_gdal
84 #else
85 INTEGER,PARAMETER :: grid_id_default = grid_id_no_driver
86 #endif
87 
88 CHARACTER(len=12),PARAMETER :: driverlist(0:2) = &
89  (/'no_driver ','grib_api ','gdal '/)
90 
91 #ifdef HAVE_LIBGDAL
92 
97  DOUBLE PRECISION :: xmin=dmiss
98  DOUBLE PRECISION :: ymin=dmiss
99  DOUBLE PRECISION :: xmax=dmiss
100  DOUBLE PRECISION :: ymax=dmiss
101 END TYPE gdal_file_id_options
102 #endif
103 
106 TYPE grid_file_id
107 PRIVATE
108 #ifdef HAVE_LIBGRIBAPI
109 INTEGER :: gaid=imiss
110 #endif
111 #ifdef HAVE_LIBGDAL
112 TYPE(gdaldataseth) :: gdalid
113 INTEGER :: nlastband=0
114 TYPE(gdal_file_id_options) :: gdal_options
115 TYPE(grid_file_id),POINTER :: file_id_copy=>null()
116 #endif
117 INTEGER :: driver=grid_id_default
118 END TYPE grid_file_id
119 
120 
123 TYPE grid_id
124 PRIVATE
125 INTEGER :: nodriverid=imiss
126 #ifdef HAVE_LIBGRIBAPI
127 INTEGER :: gaid=imiss
128 #endif
129 #ifdef HAVE_LIBGDAL
130 TYPE(gdalrasterbandh) :: gdalid
131 TYPE(grid_file_id),POINTER :: file_id=>null()
132 #endif
133 INTEGER :: driver=grid_id_default
134 END TYPE grid_id
135 
138 INTERFACE init
139  MODULE PROCEDURE grid_file_id_init, grid_id_init
140 END INTERFACE
141 
143 INTERFACE delete
144  MODULE PROCEDURE grid_file_id_delete, grid_id_delete
145 END INTERFACE
146 
148 INTERFACE copy
149  MODULE PROCEDURE grid_id_copy
150 END INTERFACE
151 
153 INTERFACE export
154  MODULE PROCEDURE grid_id_export
155 END 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 
170 
181 INTERFACE 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
183 END INTERFACE
184 
187 INTERFACE display
188  MODULE PROCEDURE grid_id_display
189 END INTERFACE
190 
191 PRIVATE 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 
194 CONTAINS
195 
196 
197 SUBROUTINE grid_file_id_init(this, filename, mode, driver, from_grid_id)
198 TYPE(grid_file_id),INTENT(out) :: this ! object to initialise
199 CHARACTER(len=*),INTENT(in) :: filename ! name of file containing gridded data, in the format [driver:]pathname
200 CHARACTER(len=*),INTENT(in) :: mode ! access mode for file, 'r' or 'w'
201 INTEGER,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
202 TYPE(grid_id),INTENT(in),OPTIONAL :: from_grid_id ! select the driver as the one associated to the provided grid_id object
203 
204 this = grid_file_id_new(filename, mode, driver, from_grid_id)
205 
206 END SUBROUTINE grid_file_id_init
207 
208 
221 FUNCTION grid_file_id_new(filename, mode, driver, from_grid_id) RESULT(this)
222 CHARACTER(len=*),INTENT(in) :: filename
223 CHARACTER(len=*),INTENT(in) :: mode
224 INTEGER,INTENT(in),OPTIONAL :: driver
225 TYPE(grid_id),INTENT(in),OPTIONAL :: from_grid_id
226 TYPE(grid_file_id) :: this
227 
228 INTEGER :: n, ier, nf
229 #ifdef HAVE_LIBGDAL
230 INTEGER :: imode
231 #endif
232 TYPE(csv_record) :: driveropts
233 CHARACTER(len=12) :: drivername
234 
235 #ifdef HAVE_LIBGDAL
236 CALL gdalnullify(this%gdalid)
237 #endif
238 
239 IF (filename == '' .OR. .NOT.c_e(filename)) RETURN
240 
241 n = index(filename,':')
242 IF (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)
271 ENDIF
272 IF (PRESENT(driver)) THEN ! override with driver
273  this%driver = driver
274 ENDIF
275 IF (PRESENT(from_grid_id)) THEN ! override with from_grid_id
276  this%driver = from_grid_id%driver
277 ENDIF
278 
279 #ifdef HAVE_LIBGRIBAPI
280 IF (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
283 ENDIF
284 #endif
285 #ifdef HAVE_LIBGDAL
286 IF (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
297 ENDIF
298 #endif
299 
300 END FUNCTION grid_file_id_new
301 
302 
306 FUNCTION grid_file_id_count(this) RESULT(count)
307 TYPE(grid_file_id),INTENT(in) :: this
308 INTEGER :: count
309 
310 INTEGER :: ier
311 
312 count = 0
313 #ifdef HAVE_LIBGRIBAPI
314 IF (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
319 ENDIF
320 #endif
321 #ifdef HAVE_LIBGDAL
322 IF (this%driver == grid_id_gdal) THEN
323  IF (gdalassociated(this%gdalid)) THEN
324  count = gdalgetrastercount(this%gdalid)
325  ENDIF
326 ENDIF
327 #endif
328 
329 END FUNCTION grid_file_id_count
330 
331 
337 SUBROUTINE grid_file_id_delete(this)
338 TYPE(grid_file_id),INTENT(inout) :: this
339 
340 #ifdef HAVE_LIBGRIBAPI
341 IF (this%driver == grid_id_grib_api) THEN
342  IF (c_e(this%gaid)) CALL grib_close_file(this%gaid)
343 ENDIF
344 this%gaid = imiss
345 #endif
346 #ifdef HAVE_LIBGDAL
347 IF (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
351 ENDIF
352 CALL gdalnullify(this%gdalid)
353 #endif
354 
355 this%driver = imiss
356 
357 END 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.
363 FUNCTION grid_file_id_c_e(this)
364 TYPE(grid_file_id),INTENT(in) :: this ! object to be checked
365 LOGICAL :: grid_file_id_c_e
366 
367 grid_file_id_c_e = .false.
368 
369 #ifdef HAVE_LIBGRIBAPI
370 IF (this%driver == grid_id_grib_api) THEN
371  grid_file_id_c_e = c_e(this%gaid)
372 ENDIF
373 #endif
374 #ifdef HAVE_LIBGDAL
375 IF (this%driver == grid_id_gdal) THEN
376  grid_file_id_c_e = gdalassociated(this%gdalid)
377 ENDIF
378 #endif
379 
380 END 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.
386 FUNCTION grid_file_id_c_e_v(this)
387 TYPE(grid_file_id),INTENT(in) :: this(:) ! object to be checked
388 LOGICAL :: grid_file_id_c_e_v(SIZE(this))
389 
390 INTEGER :: i
391 
392 DO i = 1, SIZE(this)
393  grid_file_id_c_e_v(i) = c_e(this(i))
394 ENDDO
395 
396 END FUNCTION grid_file_id_c_e_v
397 
398 
399 SUBROUTINE grid_id_init(this, from_grid_file_id, grib_api_template, grib_api_id)
400 TYPE(grid_id),INTENT(out) :: this ! object to be initialized
401 TYPE(grid_file_id),INTENT(inout),OPTIONAL :: from_grid_file_id ! file object from which grid object has to be created
402 CHARACTER(len=*),INTENT(in),OPTIONAL :: grib_api_template ! grib_api template file from which grid_object has to be created
403 INTEGER,INTENT(in),OPTIONAL :: grib_api_id ! grib_api id obtained directly from a \a grib_get subroutine call
404 
405 this = grid_id_new(from_grid_file_id, grib_api_template, grib_api_id)
406 
407 END SUBROUTINE grid_id_init
408 
409 
419 FUNCTION grid_id_new(from_grid_file_id, grib_api_template, grib_api_id, &
420  no_driver_id) RESULT(this)
421 TYPE(grid_file_id),INTENT(inout),OPTIONAL,TARGET :: from_grid_file_id
422 CHARACTER(len=*),INTENT(in),OPTIONAL :: grib_api_template
423 INTEGER,INTENT(in),OPTIONAL :: grib_api_id
424 INTEGER,INTENT(in),OPTIONAL :: no_driver_id
425 TYPE(grid_id) :: this
426 
427 INTEGER :: ier
428 
429 #ifdef HAVE_LIBGDAL
430 CALL gdalnullify(this%gdalid)
431 #endif
432 
433 IF (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
461 ELSE 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
465 ELSE IF (PRESENT(grib_api_id)) THEN
466  this%driver = grid_id_grib_api
467  this%gaid = grib_api_id
468 #endif
469 ELSE IF (PRESENT(no_driver_id)) THEN
470  this%driver = grid_id_no_driver
471  this%nodriverid = no_driver_id
472 ENDIF
473 
474 END FUNCTION grid_id_new
475 
476 
481 SUBROUTINE grid_id_delete(this)
482 TYPE(grid_id),INTENT(inout) :: this
483 
484 this%nodriverid = imiss
485 #ifdef HAVE_LIBGRIBAPI
486 IF (this%driver == grid_id_grib_api) THEN
487  IF (c_e(this%gaid)) CALL grib_release(this%gaid)
488 ENDIF
489 this%gaid = imiss
490 #endif
491 #ifdef HAVE_LIBGDAL
492 CALL gdalnullify(this%gdalid)
493 NULLIFY(this%file_id)
494 #endif
495 
496 this%driver = imiss
497 
498 END SUBROUTINE grid_id_delete
499 
500 
503 FUNCTION grid_id_readonly(this) RESULT(readonly)
504 TYPE(grid_id),INTENT(in) :: this
505 LOGICAL :: readonly
506 
507 readonly = this%driver /= grid_id_grib_api
508 
509 END FUNCTION grid_id_readonly
510 
511 
517 SUBROUTINE grid_id_copy(this, that)
518 TYPE(grid_id),INTENT(in) :: this
519 TYPE(grid_id),INTENT(out) :: that
520 
521 that = this ! start with a shallow copy
522 
523 #ifdef HAVE_LIBGRIBAPI
524 IF (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
529 ENDIF
530 #endif
531 #ifdef HAVE_LIBGDAL
532 IF (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
538 ENDIF
539 #endif
540 
541 END SUBROUTINE grid_id_copy
542 
543 
547 SUBROUTINE grid_id_export(this, file_id)
548 TYPE(grid_id),INTENT(inout) :: this
549 TYPE(grid_file_id),INTENT(in) :: file_id
550 
551 INTEGER :: ier
552 
553 IF (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
558 ENDIF
559 #ifdef HAVE_LIBGDAL
560 IF (this%driver == grid_id_gdal .AND. file_id%driver == grid_id_gdal) THEN
561  ! not implemented, log?
562 ENDIF
563 #endif
564 
565 END 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.
572 FUNCTION grid_id_c_e(this)
573 TYPE(grid_id),INTENT(in) :: this ! object to be checked
574 LOGICAL :: grid_id_c_e
575 
576 grid_id_c_e = .false.
577 
578 #ifdef HAVE_LIBGRIBAPI
579 IF (this%driver == grid_id_grib_api) THEN
580  grid_id_c_e = c_e(this%gaid)
581 ENDIF
582 #endif
583 #ifdef HAVE_LIBGDAL
584 IF (this%driver == grid_id_gdal) THEN
585  grid_id_c_e = gdalassociated(this%gdalid)
586 ENDIF
587 #endif
588 IF (this%driver == grid_id_no_driver) THEN
589  grid_id_c_e = c_e(this%nodriverid)
590 ENDIF
591 
592 END 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.
599 FUNCTION grid_id_c_e_v(this)
600 TYPE(grid_id),INTENT(in) :: this(:) ! object to be checked
601 LOGICAL :: grid_id_c_e_v(SIZE(this))
602 
603 INTEGER :: i
604 
605 DO i = 1, SIZE(this)
606  grid_id_c_e_v(i) = c_e(this(i))
607 ENDDO
608 
609 END FUNCTION grid_id_c_e_v
610 
611 
616 FUNCTION grid_file_id_get_driver(this) RESULT(driver)
617 TYPE(grid_file_id),INTENT(in) :: this
618 CHARACTER(len=LEN(driverlist)) :: driver
619 
620 IF (this%driver > 0 .AND. this%driver <= SIZE(driverlist)) THEN
621  driver = driverlist(this%driver)
622 ELSE
623  driver = driverlist(0)
624 ENDIF
625 
626 END FUNCTION grid_file_id_get_driver
627 
628 
633 FUNCTION grid_id_get_driver(this) RESULT(driver)
634 TYPE(grid_id),INTENT(in) :: this
635 CHARACTER(len=LEN(driverlist)) :: driver
636 
637 IF (this%driver > 0 .AND. this%driver <= SIZE(driverlist)) THEN
638  driver = driverlist(this%driver)
639 ELSE
640  driver = driverlist(0)
641 ENDIF
642 
643 END FUNCTION grid_id_get_driver
644 
645 
652 SUBROUTINE grid_id_display(this, namespace)
653 TYPE(grid_id),INTENT(in) :: this
654 CHARACTER(len=*),OPTIONAL :: namespace
655 
656 INTEGER :: kiter, iret
657 CHARACTER(len=255) :: key, value, lnamespace
658 
659 
660 #ifdef HAVE_LIBGRIBAPI
661 IF (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)
669 
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 
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 
693 ENDIF
694 
695 #endif
696 
697 END SUBROUTINE grid_id_display
698 
699 
700 #ifdef HAVE_LIBGRIBAPI
701 
703 FUNCTION grid_file_id_get_gaid(this) RESULT(gaid)
704 TYPE(grid_file_id),INTENT(in) :: this
705 INTEGER :: gaid
706 gaid = this%gaid
707 END FUNCTION grid_file_id_get_gaid
708 
711 FUNCTION grid_id_get_gaid(this) RESULT(gaid)
712 TYPE(grid_id),INTENT(in) :: this
713 INTEGER :: gaid
714 gaid = this%gaid
715 END FUNCTION grid_id_get_gaid
716 #endif
717 
718 
719 #ifdef HAVE_LIBGDAL
720 
722 FUNCTION grid_file_id_get_gdalid(this) RESULT(gdalid)
723 TYPE(grid_file_id),INTENT(in) :: this
724 TYPE(gdaldataseth) :: gdalid
725 gdalid = this%gdalid
726 END FUNCTION grid_file_id_get_gdalid
727 
730 FUNCTION grid_id_get_gdalid(this) RESULT(gdalid)
731 TYPE(grid_id),INTENT(in) :: this
732 TYPE(gdalrasterbandh) :: gdalid
733 gdalid = this%gdalid
734 END FUNCTION grid_id_get_gdalid
735 
738 FUNCTION grid_id_get_gdal_options(this) RESULT(gdal_options)
739 TYPE(grid_id),INTENT(in) :: this
740 TYPE(gdal_file_id_options) :: gdal_options
741 
742 TYPE(gdal_file_id_options) :: gdal_options_local
743 
744 IF (ASSOCIATED(this%file_id)) THEN
745  gdal_options = this%file_id%gdal_options
746 ELSE
747  gdal_options = gdal_options_local ! empty object
748 ENDIF
749 
750 END FUNCTION grid_id_get_gdal_options
751 #endif
752 
753 
757 SUBROUTINE grid_id_decode_data(this, field)
758 TYPE(grid_id),INTENT(in) :: this
759 REAL,INTENT(out) :: field(:,:)
760 
761 LOGICAL :: done
762 
763 done = .false.
764 #ifdef HAVE_LIBGRIBAPI
765 IF (c_e(this%gaid)) THEN
766  CALL grid_id_decode_data_gribapi(this%gaid, field)
767  done = .true.
768 ENDIF
769 #endif
770 #ifdef HAVE_LIBGDAL
771 ! subarea?
772 IF (gdalassociated(this%gdalid)) THEN
773  CALL grid_id_decode_data_gdal(this%gdalid, field, this%file_id%gdal_options)
774  done = .true.
775 ENDIF
776 #endif
777 IF (.NOT.done) field(:,:) = rmiss
778 
779 END SUBROUTINE grid_id_decode_data
780 
781 
785 SUBROUTINE grid_id_encode_data(this, field)
786 TYPE(grid_id),INTENT(inout) :: this
787 REAL,intent(in) :: field(:,:)
788 
789 #ifdef HAVE_LIBGRIBAPI
790 IF (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)
797 ENDIF
798 #endif
799 #ifdef HAVE_LIBGDAL
800 IF (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?
804 ENDIF
805 #endif
806 
807 END SUBROUTINE grid_id_encode_data
808 
809 
810 #ifdef HAVE_LIBGRIBAPI
811 SUBROUTINE grid_id_decode_data_gribapi(gaid, field)
812 INTEGER,INTENT(in) :: gaid ! grib_api id
813 REAL,INTENT(out) :: field(:,:) ! array of decoded values
814 
815 INTEGER :: editionnumber
816 INTEGER :: alternativerowscanning, &
817  iscansnegatively, jscanspositively, jpointsareconsecutive
818 INTEGER :: numberofvalues,numberofpoints
819 REAL :: vector(size(field))
820 INTEGER :: x1, x2, xs, y1, y2, ys, ord(2), ierr
821 
822 
823 call grib_get(gaid,'GRIBEditionNumber',editionnumber)
824 
825 if (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 
836 else 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 
844 end if
845 
846 CALL grib_get(gaid,'iScansNegatively',iscansnegatively,ierr)
847 IF (ierr /= grib_success) iscansnegatively=0
848 CALL grib_get(gaid,'jScansPositively',jscanspositively,ierr)
849 IF (ierr /= grib_success) jscanspositively=1
850 CALL grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive,ierr)
851 IF (ierr /= grib_success) jpointsareconsecutive=0
852 
853 call grib_get(gaid,'numberOfPoints',numberofpoints)
854 call grib_get(gaid,'numberOfValues',numberofvalues)
855 
856 IF (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
864 ENDIF
865 
866 ! get data values
867 #ifdef DEBUG
868 call l4f_log(l4f_info,'grib_api number of values: '//to_char(numberofvalues))
869 call l4f_log(l4f_info,'grib_api number of points: '//to_char(numberofpoints))
870 #endif
871 
872 CALL grib_set(gaid,'missingValue',rmiss)
873 CALL grib_get(gaid,'values',vector)
874 ! suspect bug in grib_api, when all field is missing it is set to zero
875 IF (numberofvalues == 0) vector = rmiss
876 
877 #ifdef DEBUG
878 CALL 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))))
880 CALL l4f_log(l4f_debug, 'grib_api, decoded field with number of missing: '// &
881  t2c(count(.NOT.c_e(vector))))
882 #endif
883 
884 IF (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()
889 ENDIF
890 
891 ! Transfer data field changing scanning mode to 64
892 IF (iscansnegatively == 0) THEN
893  x1 = 1
894  x2 = SIZE(field,1)
895  xs = 1
896 ELSE
897  x1 = SIZE(field,1)
898  x2 = 1
899  xs = -1
900 ENDIF
901 IF (jscanspositively == 0) THEN
902  y1 = SIZE(field,2)
903  y2 = 1
904  ys = -1
905 ELSE
906  y1 = 1
907  y2 = SIZE(field,2)
908  ys = 1
909 ENDIF
910 
911 IF ( jpointsareconsecutive == 0) THEN
912  ord = (/1,2/)
913 ELSE
914  ord = (/2,1/)
915 ENDIF
916 
917 field(x1:x2:xs,y1:y2:ys) = reshape(vector, &
918  (/SIZE(field,1),SIZE(field,2)/), order=ord)
919 
920 END SUBROUTINE grid_id_decode_data_gribapi
921 
922 
923 SUBROUTINE grid_id_encode_data_gribapi(gaid, field)
924 INTEGER,INTENT(in) :: gaid ! grib_api id
925 REAL,intent(in) :: field(:,:) ! data array to be encoded
926 
927 INTEGER :: editionnumber
928 INTEGER :: alternativerowscanning, iscansnegatively, &
929  jscanspositively, jpointsareconsecutive
930 INTEGER :: x1, x2, xs, y1, y2, ys, ierr
931 
932 call grib_get(gaid,'GRIBEditionNumber',editionnumber)
933 
934 if (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 
944 else if( editionnumber /= 1) then
945 
946  call l4f_log(l4f_error, &
947  "grib_api GribEditionNumber not supported: "//t2c(editionnumber))
948  call raise_error()
949  RETURN
950 
951 end if
952 
953 CALL grib_get(gaid,'iScansNegatively',iscansnegatively,ierr)
954 IF (ierr /= grib_success) iscansnegatively=0
955 CALL grib_get(gaid,'jScansPositively',jscanspositively,ierr)
956 IF (ierr /= grib_success) jscanspositively=1
957 CALL grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive,ierr)
958 IF (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
964 CALL l4f_log(l4f_debug, 'grib_api, Ni,Nj:'//t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
965 #endif
966 CALL grib_set(gaid,'Ni',SIZE(field,1), ierr)
967 CALL grib_set(gaid,'Nj',SIZE(field,2), ierr)
968 
969 ! Transfer data field changing scanning mode from 64
970 IF (iscansnegatively == 0) THEN
971  x1 = 1
972  x2 = SIZE(field,1)
973  xs = 1
974 ELSE
975  x1 = SIZE(field,1)
976  x2 = 1
977  xs = -1
978 ENDIF
979 IF (jscanspositively == 0) THEN
980  y1 = SIZE(field,2)
981  y2 = 1
982  ys = -1
983 ELSE
984  y1 = 1
985  y2 = SIZE(field,2)
986  ys = 1
987 ENDIF
988 
989 
990 IF (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 
1004 ELSE
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 
1014 ENDIF
1015 
1016 #ifdef DEBUG
1017 CALL 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))))
1019 CALL l4f_log(l4f_debug, 'grib_api, coding field with number of missing: '// &
1020  t2c(count(.NOT.c_e(field))))
1021 CALL l4f_log(l4f_debug, 'grib_api, sizex:'//t2c(x1)//','//t2c(x2)//','//t2c(xs))
1022 CALL l4f_log(l4f_debug, 'grib_api, sizey:'//t2c(y1)//','//t2c(y2)//','//t2c(ys))
1023 #endif
1024 IF (jpointsareconsecutive == 0) THEN
1025  CALL grib_set(gaid,'values', reshape(field(x1:x2:xs,y1:y2:ys), &
1026  (/SIZE(field)/)))
1027 ELSE
1028  CALL grib_set(gaid,'values', reshape(transpose(field(x1:x2:xs,y1:y2:ys)), &
1029  (/SIZE(field)/)))
1030 ENDIF
1031 
1032 END SUBROUTINE grid_id_encode_data_gribapi
1033 #endif
1034 
1035 
1036 #ifdef HAVE_LIBGDAL
1037 SUBROUTINE grid_id_decode_data_gdal(gdalid, field, gdal_options)
1038 #ifdef F2003_FULL_FEATURES
1039 USE ieee_arithmetic
1040 #endif
1041 TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal id
1042 REAL,INTENT(out) :: field(:,:) ! array of decoded values
1043 TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
1044 
1045 TYPE(gdaldataseth) :: hds
1046 REAL(kind=c_double) :: geotrans(6), dummy1, dummy2, dummy3, dummy4
1047 REAL :: gdalmiss
1048 REAL,ALLOCATABLE :: buffer(:,:)
1049 INTEGER :: ix1, iy1, ix2, iy2, ixs, iys, ord(2), ier
1050 INTEGER(kind=c_int) :: nrx, nry
1051 LOGICAL :: must_trans
1052 
1053 
1054 hds = gdalgetbanddataset(gdalid) ! go back to dataset
1055 ier = gdalgetgeotransform(hds, geotrans)
1056 
1057 IF (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 
1083 ELSE 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 
1109 ELSE ! 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
1114 ENDIF
1115 
1116 ! read data from file
1117 CALL gdalrastersimpleread_f(gdalid, gdal_options%xmin, gdal_options%ymin, &
1118  gdal_options%xmax, gdal_options%ymax, buffer, dummy1, dummy2, dummy3, dummy4)
1119 
1120 IF (.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
1125 ENDIF
1126 
1127 IF (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
1135 ENDIF
1136 
1137 #ifdef F2003_FULL_FEATURES
1138 ! aparently gdal datasets may contain NaN
1139 WHERE(ieee_is_nan(buffer))
1140  buffer = rmiss
1141 END WHERE
1142 #else
1143 WHERE(buffer /= buffer)
1144  buffer = rmiss
1145 END WHERE
1146 #endif
1147 
1148 ! set missing value if necessary
1149 gdalmiss = real(gdalgetrasternodatavalue(gdalid, ier))
1150 IF (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
1157 ELSE
1158 #ifdef DEBUG
1159  CALL l4f_log(l4f_info, 'gdal no missing data found in band')
1160 #endif
1161 ENDIF
1162 
1163 ! reshape the field
1164 IF (must_trans) THEN
1165  field(ix1:ix2:ixs,iy1:iy2:iys) = transpose(buffer)
1166 ELSE
1167  field(ix1:ix2:ixs,iy1:iy2:iys) = buffer(:,:)
1168 ENDIF
1169 
1170 
1171 END SUBROUTINE grid_id_decode_data_gdal
1172 #endif
1173 
1174 END 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.
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.