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
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
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 
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
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)
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
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.