libsim  Versione7.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
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
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)
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
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
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
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
Definitions of constants and functions for working with missing values.
Functions that return a trimmed CHARACTER representation of the input variable.
Gestione degli errori.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
This module defines an abstract interface to different drivers for access to files containing gridded...
Derived type containing driver-specific options for gdal.
Derived type associated to a file-like object containing many blocks/messages/records/bands of gridde...
Destructors for the corresponding classes.
Methods for successively obtaining the fields of a csv_record object.
Check whether the corresponding object has been correctly associated.
Index method.
Restituiscono il valore dell&#39;oggetto in forma di stringa stampabile.
Make a deep copy, if possible, of the grid identifier.
Constructors for the corresponding classes in SUBROUTINE form.
Class for interpreting the records of a csv file.
Export a grid to a file.
classe per la gestione del logging
Utilities for CHARACTER variables.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Utilities for managing files.
Display on standard output a description of the grid_id object provided.

Generated with Doxygen.