libsim  Versione6.3.0
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
310 INTEGER :: ier
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
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)
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
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
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
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)
821 
822 
823 call grib_get(gaid,'GRIBEditionNumber',editionnumber)
824 
825 if (editionnumber == 2) then
826 
827  call grib_get(gaid,'alternativeRowScanning',alternativerowscanning)
828  if (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  end if
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)
847 call grib_get(gaid,'jScansPositively',jscanspositively)
848 call grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive)
849 
850 call grib_get(gaid,'numberOfPoints',numberofpoints)
851 call grib_get(gaid,'numberOfValues',numberofvalues)
852 
853 IF (numberofpoints /= SIZE(field)) THEN
854  CALL l4f_log(l4f_error, 'grid_id_decode_data_gribapi numberOfPoints and grid size different')
855  CALL l4f_log(l4f_error, 'grid_id_decode_data_gribapi numberOfPoints: ' &
856  //t2c(numberofpoints)//', nx,ny: '&
857  //t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
858  CALL raise_error()
859  field(:,:) = rmiss
860  RETURN
861 ENDIF
862 
863 ! get data values
864 #ifdef DEBUG
865 call l4f_log(l4f_info,'grib_api number of values: '//to_char(numberofvalues))
866 call l4f_log(l4f_info,'grib_api number of points: '//to_char(numberofpoints))
867 #endif
868 
869 CALL grib_set(gaid,'missingValue',rmiss)
870 CALL grib_get(gaid,'values',vector)
871 ! suspect bug in grib_api, when all field is missing it is set to zero
872 IF (numberofvalues == 0) vector = rmiss
873 
874 #ifdef DEBUG
875 CALL l4f_log(l4f_debug, 'grib_api, decoded field in interval: '// &
876  t2c(minval(vector,mask=c_e(vector)))//' '//t2c(maxval(vector,mask=c_e(vector))))
877 CALL l4f_log(l4f_debug, 'grib_api, decoded field with number of missing: '// &
878  t2c(count(.NOT.c_e(vector))))
879 #endif
880 
881 IF (numberofvalues /= count(c_e(vector))) THEN
882  CALL l4f_log(l4f_warn, 'grid_id_decode_data_gribapi numberOfValues and valid data count different')
883  CALL l4f_log(l4f_warn, 'grid_id_decode_data_gribapi numberOfValues: ' &
884  //t2c(numberofvalues)//', valid data: '//t2c(count(c_e(vector))))
885 ! CALL raise_warning()
886 ENDIF
887 
888 ! Transfer data field changing scanning mode to 64
889 IF (iscansnegatively == 0) THEN
890  x1 = 1
891  x2 = SIZE(field,1)
892  xs = 1
893 ELSE
894  x1 = SIZE(field,1)
895  x2 = 1
896  xs = -1
897 ENDIF
898 IF (jscanspositively == 0) THEN
899  y1 = SIZE(field,2)
900  y2 = 1
901  ys = -1
902 ELSE
903  y1 = 1
904  y2 = SIZE(field,2)
905  ys = 1
906 ENDIF
907 
908 IF ( jpointsareconsecutive == 0) THEN
909  ord = (/1,2/)
910 ELSE
911  ord = (/2,1/)
912 ENDIF
913 
914 field(x1:x2:xs,y1:y2:ys) = reshape(vector, &
915  (/SIZE(field,1),SIZE(field,2)/), order=ord)
916 
917 END SUBROUTINE grid_id_decode_data_gribapi
918 
919 
920 SUBROUTINE grid_id_encode_data_gribapi(gaid, field)
921 INTEGER,INTENT(in) :: gaid ! grib_api id
922 REAL,intent(in) :: field(:,:) ! data array to be encoded
923 
924 INTEGER :: editionnumber
925 INTEGER :: alternativerowscanning, iscansnegatively, &
926  jscanspositively, jpointsareconsecutive
927 INTEGER :: x1, x2, xs, y1, y2, ys
928 
929 
930 call grib_get(gaid,'GRIBEditionNumber',editionnumber)
931 
932 if (editionnumber == 2) then
933 
934  call grib_get(gaid,'alternativeRowScanning',alternativerowscanning)
935  if (alternativerowscanning /= 0)then
936  call l4f_log(l4f_error, "grib_api alternativeRowScanning not supported: " &
937  //trim(to_char(alternativerowscanning)))
938  call raise_error()
939  RETURN
940  end if
941 
942 else if( editionnumber /= 1) then
943 
944  call l4f_log(l4f_error, &
945  "grib_api GribEditionNumber not supported: "//t2c(editionnumber))
946  call raise_error()
947  RETURN
948 
949 end if
950 
951 call grib_get(gaid,'iScansNegatively',iscansnegatively)
952 call grib_get(gaid,'jScansPositively',jscanspositively)
953 call grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive)
954 
955 ! queste sono gia` fatte in export_gridinfo, si potrebbero evitare?!
956 #ifdef DEBUG
957 CALL l4f_log(l4f_debug, 'grib_api, Ni,Nj:'//t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
958 #endif
959 call grib_set(gaid,'Ni',SIZE(field,1))
960 call grib_set(gaid,'Nj',SIZE(field,2))
961 
962 ! Transfer data field changing scanning mode from 64
963 IF (iscansnegatively == 0) THEN
964  x1 = 1
965  x2 = SIZE(field,1)
966  xs = 1
967 ELSE
968  x1 = SIZE(field,1)
969  x2 = 1
970  xs = -1
971 ENDIF
972 IF (jscanspositively == 0) THEN
973  y1 = SIZE(field,2)
974  y2 = 1
975  ys = -1
976 ELSE
977  y1 = 1
978  y2 = SIZE(field,2)
979  ys = 1
980 ENDIF
981 
982 
983 if (any(field == rmiss)) then
984 
985  call grib_set(gaid,'missingValue',rmiss)
986  if (editionnumber == 1) then
987 ! enable bitmap in grib1
988 ! grib_api 1.9.9 goes into an infinite loop with second order packing here
989  CALL grib_set(gaid,'packingType','grid_simple')
990  call grib_set(gaid,"bitmapPresent",1)
991  else
992 ! enable bitmap in grib2
993  call grib_set(gaid,"bitMapIndicator",0)
994  endif
995 
996 else
997 
998  if (editionnumber == 1) then
999 ! disable bitmap in grib1
1000  call grib_set(gaid,"bitmapPresent",0)
1001  else
1002 ! disable bitmap in grib2
1003  call grib_set(gaid,"bitMapIndicator",255)
1004  endif
1005 
1006 end if
1007 
1008 !TODO: gestire il caso TUTTI dati mancanti
1009 
1010 #ifdef DEBUG
1011 CALL l4f_log(l4f_debug, 'grib_api, coding field in interval: '// &
1012  t2c(minval(field,mask=c_e(field)))//' '//t2c(maxval(field,mask=c_e(field))))
1013 CALL l4f_log(l4f_debug, 'grib_api, coding field with number of missing: '// &
1014  t2c(count(.NOT.c_e(field))))
1015 CALL l4f_log(l4f_debug, 'grib_api, sizex:'//t2c(x1)//','//t2c(x2)//','//t2c(xs))
1016 CALL l4f_log(l4f_debug, 'grib_api, sizey:'//t2c(y1)//','//t2c(y2)//','//t2c(ys))
1017 #endif
1018 IF ( jpointsareconsecutive == 0) THEN
1019  CALL grib_set(gaid,'values', reshape(field(x1:x2:xs,y1:y2:ys), &
1020  (/SIZE(field)/)))
1021 ELSE
1022  CALL grib_set(gaid,'values', reshape(transpose(field(x1:x2:xs,y1:y2:ys)), &
1023  (/SIZE(field)/)))
1024 ENDIF
1025 
1026 END SUBROUTINE grid_id_encode_data_gribapi
1027 #endif
1028 
1029 
1030 #ifdef HAVE_LIBGDAL
1031 SUBROUTINE grid_id_decode_data_gdal(gdalid, field, gdal_options)
1032 TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal id
1033 REAL,INTENT(out) :: field(:,:) ! array of decoded values
1034 TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
1035 
1036 TYPE(gdaldataseth) :: hds
1037 REAL(kind=c_double) :: geotrans(6), dummy1, dummy2, dummy3, dummy4
1038 REAL :: gdalmiss
1039 REAL,ALLOCATABLE :: buffer(:,:)
1040 INTEGER :: ix1, iy1, ix2, iy2, ixs, iys, ord(2), ier
1041 INTEGER(kind=c_int) :: nrx, nry
1042 LOGICAL :: must_trans
1043 
1044 
1045 hds = gdalgetbanddataset(gdalid) ! go back to dataset
1046 ier = gdalgetgeotransform(hds, geotrans)
1047 
1048 IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN
1049 ! transformation is diagonal, no transposing
1050  IF (geotrans(2) > 0.0_c_double) THEN
1051  ix1 = 1
1052  ix2 = SIZE(field,1)
1053  ixs = 1
1054  ELSE
1055  ix1 = SIZE(field,1)
1056  ix2 = 1
1057  ixs = -1
1058  ENDIF
1059  IF (geotrans(6) > 0.0_c_double) THEN
1060  iy1 = 1
1061  iy2 = SIZE(field,2)
1062  iys = 1
1063  ELSE
1064  iy1 = SIZE(field,2)
1065  iy2 = 1
1066  iys = -1
1067  ENDIF
1068  nrx = SIZE(field,1)
1069  nry = SIZE(field,2)
1070  ord = (/1,2/)
1071  must_trans = .false.
1072 ! ALLOCATE(buffer(nrx,nry))
1073 
1074 ELSE IF (geotrans(2) == 0.0_c_double .AND. geotrans(6) == 0.0_c_double) THEN
1075 ! transformation is anti-diagonal, transposing required this should not happen
1076  IF (geotrans(3) > 0.0_c_double) THEN
1077  ix1 = 1
1078  ix2 = SIZE(field,1)
1079  ixs = 1
1080  ELSE
1081  ix1 = SIZE(field,1)
1082  ix2 = 1
1083  ixs = -1
1084  ENDIF
1085  IF (geotrans(5) > 0.0_c_double) THEN
1086  iy1 = 1
1087  iy2 = SIZE(field,2)
1088  iys = 1
1089  ELSE
1090  iy1 = SIZE(field,2)
1091  iy2 = 1
1092  iys = -1
1093  ENDIF
1094  nrx = SIZE(field,2)
1095  nry = SIZE(field,1)
1096  ord = (/2,1/)
1097  must_trans = .true.
1098 ! ALLOCATE(buffer(nry,nrx))
1099 
1100 ELSE ! transformation is a rotation, not supported
1101  CALL l4f_log(l4f_error, 'gdal geotransform is a generic rotation, not supported')
1102  CALL raise_error()
1103  field(:,:) = rmiss
1104  RETURN
1105 ENDIF
1106 
1107 ! read data from file
1108 CALL gdalrastersimpleread_f(gdalid, gdal_options%xmin, gdal_options%ymin, &
1109  gdal_options%xmax, gdal_options%ymax, buffer, dummy1, dummy2, dummy3, dummy4)
1110 
1111 IF (.NOT.ALLOCATED(buffer)) THEN ! error in read
1112  CALL l4f_log(l4f_error, 'gdal error in reading with gdal driver')
1113  CALL raise_error()
1114  field(:,:) = rmiss
1115  RETURN
1116 ENDIF
1117 
1118 IF (SIZE(buffer) /= (SIZE(field)))THEN
1119  CALL l4f_log(l4f_error, 'gdal raster band and gridinfo size different')
1120  CALL l4f_log(l4f_error, 'gdal rasterband: ' &
1121  //t2c(SIZE(buffer,1))//'X'//t2c(SIZE(buffer,2))//', nx,ny:' &
1122  //t2c(SIZE(field,ord(1)))//'X'//t2c(SIZE(field,ord(2))))
1123  CALL raise_error()
1124  field(:,:) = rmiss
1125  RETURN
1126 ENDIF
1127 
1128 ! set missing value if necessary
1129 gdalmiss = REAL(gdalgetrasternodatavalue(gdalid, ier))
1130 IF (ier /= 0) THEN ! success -> there are missing values
1131 #ifdef DEBUG
1132  CALL l4f_log(l4f_info, 'gdal missing data value: '//trim(to_char(gdalmiss)))
1133 #endif
1134  WHERE(buffer(:,:) == gdalmiss)
1135  buffer(:,:) = rmiss
1136  END WHERE
1137 ELSE
1138 #ifdef DEBUG
1139  CALL l4f_log(l4f_info, 'gdal no missing data found in band')
1140 #endif
1141 ENDIF
1142 
1143 ! reshape the field
1144 IF (must_trans) THEN
1145  field(ix1:ix2:ixs,iy1:iy2:iys) = transpose(buffer)
1146 ELSE
1147  field(ix1:ix2:ixs,iy1:iy2:iys) = buffer(:,:)
1148 ENDIF
1149 
1150 
1151 END SUBROUTINE grid_id_decode_data_gdal
1152 #endif
1153 
1154 END MODULE grid_id_class
Set of functions that return a trimmed CHARACTER representation of the input variable.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
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...
This module defines an abstract interface to different drivers for access to files containing gridded...
Destructors for the corresponding classes.
Utilities for managing files.
Methods for successively obtaining the fields of a csv_record object.
Check whether the corresponding object has been correctly associated.
Index method.
Set of functions that return a CHARACTER representation of the input variable.
Gestione degli errori.
Make a deep copy, if possible, of the grid identifier.
Definitions of constants and functions for working with missing values.
Constructors for the corresponding classes in SUBROUTINE form.
Export a grid to a file.
classe per la gestione del logging
Utilities for CHARACTER variables.
Display on standard output a description of the grid_id object provided.

Generated with Doxygen.