gdal.F90

!    Copyright 2011 Davide Cesari <dcesari69 at gmail dot com>
!
!    This file is part of FortranGIS.
!
!    FortranGIS is free software: you can redistribute it and/or modify
!    it under the terms of the GNU Lesser General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    FortranGIS is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public
!    License along with FortranGIS.  If not, see
!    <http://www.gnu.org/licenses/>.

#include "gdalproto_doxy.f90"
!! 
!! As a general guideline, note that when a \c char** object is
!! encountered in the C interface, it should usually be interfaced in
!! Fortran by means of the fortranc::c_ptr_ptr derived type.
!!
!! Other Fortran-style subroutines, functions and procedure interfaces
!! are documented explicitely here.
!!
!! For an example of application of the \a gdal module, please refer
!! to the following test program, which creates a very simple gdal
!! raster dataset, exports it on a GEOTiff file and successively reads
!! it:
!! \include gdal_test.F90
!!
!! \ingroup libfortrangis
MODULE gdal
use,INTRINSIC :: iso_c_binding
IMPLICIT NONE

! Hand made symbolic constant definitions
! GDALDataType
INTEGER(kind=c_int),PARAMETER :: gdt_unknown = 0
INTEGER(kind=c_int),PARAMETER :: gdt_byte = 1
INTEGER(kind=c_int),PARAMETER :: gdt_uint16 = 2
INTEGER(kind=c_int),PARAMETER :: gdt_int16 = 3
INTEGER(kind=c_int),PARAMETER :: gdt_uint32 = 4
INTEGER(kind=c_int),PARAMETER :: gdt_int32 = 5
INTEGER(kind=c_int),PARAMETER :: gdt_float32 = 6
INTEGER(kind=c_int),PARAMETER :: gdt_float64 = 7
INTEGER(kind=c_int),PARAMETER :: gdt_cint16 = 8
INTEGER(kind=c_int),PARAMETER :: gdt_cint32 = 9
INTEGER(kind=c_int),PARAMETER :: gdt_cfloat32 = 10
INTEGER(kind=c_int),PARAMETER :: gdt_cfloat64 = 11
INTEGER(kind=c_int),PARAMETER :: gdt_typecount = 12

! GDALAccess
INTEGER(kind=c_int),PARAMETER :: ga_readonly = 0
INTEGER(kind=c_int),PARAMETER :: ga_update = 1

! GDALRWFlag
INTEGER(kind=c_int),PARAMETER :: gf_read = 0
INTEGER(kind=c_int),PARAMETER :: gf_write = 1

INTEGER(kind=c_int),PARAMETER :: & ! GDALColorInterp
 gci_undefined = 0, gci_grayindex = 1, gci_paletteindex = 2, &
 gci_redband = 3, gci_greenband = 4, gci_blueband = 5, &
 gci_alphaband = 6, gci_hueband = 7, gci_saturationband = 8, &
 gci_lightnessband = 9, gci_cyanband = 10, gci_magentaband = 11, &
 gci_yellowband = 12, gci_blackband = 13, gci_ycbcr_yband = 14, &
 gci_ycbcr_cbband = 15,gci_ycbcr_crband = 16, gci_max = 16

INTEGER(kind=c_int),PARAMETER :: & ! GDALPaletteInterp
 gpi_gray = 0, gpi_rgb = 1, gpi_cmyk = 2, gpi_hls = 3

INTEGER(kind=c_int),PARAMETER :: & ! GDALRATFieldType
 gft_integer = 0, gft_real = 1, gft_string = 2

INTEGER(kind=c_int),PARAMETER :: & ! GDALRATFieldUsage
 gfu_generic = 0, gfu_pixelcount = 1, gfu_name = 2, gfu_min = 3, &
 gfu_max = 4, gfu_minmax = 5, gfu_red = 6, gfu_green = 7, &
 gfu_blue = 8, gfu_alpha = 9, gfu_redmin = 10, gfu_greenmin = 11, &
 gfu_bluemin = 12, gfu_alphamin = 13, gfu_redmax = 14, &
 gfu_greenmax = 15, gfu_bluemax = 16, gfu_alphamax = 17, gfu_maxcount = 18

! Hand made type definitions strictly reflecting C definitions
TYPE,bind(c) :: gdal_gcp
  TYPE(c_ptr) :: pszid, pszinfo
  REAL(kind=c_double) :: dfgcppixel, dfgcpline, dfgcpx, dfgcpy, dfgcpz
END TYPE gdal_gcp

TYPE,bind(c) :: gdalrpcinfo
  REAL(kind=c_double) :: dfline_off, dfsamp_off, dflat_off, dflong_off, dfheight_off
  REAL(kind=c_double) :: dfline_scale, dfsamp_scale, dflat_scale, dflong_scale, dfheight_scale
  REAL(kind=c_double) :: adfline_num_coeff(20), adfline_den_coeff(20), &
   adfsamp_num_coeff(20), adfsamp_den_coeff(20)
  REAL(kind=c_double) :: dfmin_long, dfmin_lat, dfmax_long, dfmax_lat
END TYPE gdalrpcinfo

TYPE,bind(c) :: gdalcolorentry
  INTEGER(kind=c_short) :: c1, c2, c3, c4
END TYPE gdalcolorentry

! Machine made type definitions
include 'gdalproto_type.f90'

! Hand made interface definitions
INTERFACE
  SUBROUTINE gdalapplygeotransform(padfgeotransform, dfpixel, dfline, pdfgeox, pdfgeoy) BIND(C,name='GDALApplyGeoTransform')
  import
!GCC$ ATTRIBUTES STDCALL :: GDALApplyGeoTransform
  REAL(kind=c_double) :: padfgeotransform(*)
  REAL(kind=c_double),VALUE :: dfpixel
  REAL(kind=c_double),VALUE :: dfline
  REAL(kind=c_double) :: pdfgeox
  REAL(kind=c_double) :: pdfgeoy
  END SUBROUTINE gdalapplygeotransform
END INTERFACE

INTERFACE
  FUNCTION gdalgcpstogeotransform(ngcpcount, pasgcps, padfgeotransform, bapproxok) &
   bind(c,name='GDALGCPsToGeoTransform')
  import
!GCC$ ATTRIBUTES STDCALL :: GDALGCPsToGeoTransform
  INTEGER(kind=c_int),VALUE :: ngcpcount
  TYPE(gdal_gcp),INTENT(in) :: pasgcps
  REAL(kind=c_double) :: padfgeotransform(*)
  INTEGER(kind=c_int),VALUE :: bapproxok
  INTEGER(kind=c_int) :: gdalgcpstogeotransform
  END FUNCTION gdalgcpstogeotransform
END INTERFACE

! Machine made interface definitions
include 'gdalproto_interf.f90'

! Fortran style interfaces

INTERFACE gdalapplygeotransform_f
  MODULE PROCEDURE gdalapplygeotransform_f_0d, gdalapplygeotransform_f_1d, &
   gdalapplygeotransform_f_2d, gdalapplygeotransform_f_3d
END INTERFACE

PRIVATE gdalapplygeotransform_f_0d, gdalapplygeotransform_f_1d, &
 gdalapplygeotransform_f_2d, gdalapplygeotransform_f_3d


INTERFACE gdaldatasetrasterio_f
  MODULE PROCEDURE gdaldatasetrasterio_int8, gdaldatasetrasterio_int16, &
   gdaldatasetrasterio_int32, &
   gdaldatasetrasterio_float, gdaldatasetrasterio_double, &
   gdaldatasetrasterio_float_cmplx, gdaldatasetrasterio_double_cmplx
END INTERFACE

PRIVATE gdaldatasetrasterio_int8, gdaldatasetrasterio_int16, &
 gdaldatasetrasterio_int32, &
 gdaldatasetrasterio_float, gdaldatasetrasterio_double, &
 gdaldatasetrasterio_float_cmplx, gdaldatasetrasterio_double_cmplx

INTERFACE gdalrasterio_f
  MODULE PROCEDURE gdalrasterio_int8, gdalrasterio_int16, &
   gdalrasterio_int32, &
   gdalrasterio_float, gdalrasterio_double, &
   gdalrasterio_float_cmplx, gdalrasterio_double_cmplx
END INTERFACE

PRIVATE gdalrasterio_int8, gdalrasterio_int16, &
 gdalrasterio_int32, &
 gdalrasterio_float, gdalrasterio_double, &
 gdalrasterio_float_cmplx, gdalrasterio_double_cmplx

INTERFACE gdalmajorobjecth_new
  MODULE PROCEDURE gdalmajorobject_fromdataset_new, &
   gdalmajorobject_fromrasterband_new, &
   gdalmajorobject_fromdriver_new
END INTERFACE gdalmajorobjecth_new

PRIVATE gdalmajorobject_fromdataset_new, &
 gdalmajorobject_fromrasterband_new, &
 gdalmajorobject_fromdriver_new

! internal interfaces
INTERFACE gdaldatasetrasterio_loc
  MODULE PROCEDURE gdaldatasetrasterio_int8_loc, gdaldatasetrasterio_int16_loc, &
   gdaldatasetrasterio_int32_loc, &
   gdaldatasetrasterio_float_loc, gdaldatasetrasterio_double_loc, &
   gdaldatasetrasterio_float_cmplx_loc, gdaldatasetrasterio_double_cmplx_loc
END INTERFACE
PRIVATE gdaldatasetrasterio_loc
PRIVATE gdaldatasetrasterio_int8_loc, gdaldatasetrasterio_int16_loc, &
 gdaldatasetrasterio_int32_loc, &
 gdaldatasetrasterio_float_loc, gdaldatasetrasterio_double_loc, &
 gdaldatasetrasterio_float_cmplx_loc, gdaldatasetrasterio_double_cmplx_loc

INTERFACE gdalrasterio_loc
  MODULE PROCEDURE gdalrasterio_int8_loc, gdalrasterio_int16_loc, &
   gdalrasterio_int32_loc, &
   gdalrasterio_float_loc, gdalrasterio_double_loc, &
   gdalrasterio_float_cmplx_loc, gdalrasterio_double_cmplx_loc
END INTERFACE
PRIVATE gdalrasterio_loc
PRIVATE gdalrasterio_int8_loc, gdalrasterio_int16_loc, &
 gdalrasterio_int32_loc, &
 gdalrasterio_float_loc, gdalrasterio_double_loc, &
 gdalrasterio_float_cmplx_loc, gdalrasterio_double_cmplx_loc

CONTAINS

! Machine made procedure definitions
include 'gdalproto_proc.f90'

! Fortran specific version of some functions
FUNCTION gdalgcpstogeotransform_f(pasgcps, padfgeotransform, bapproxok)
TYPE(gdal_gcp),INTENT(in) :: pasgcps(:)
REAL(kind=c_double),INTENT(out) :: padfgeotransform(6)
INTEGER(kind=c_int),VALUE :: bapproxok
INTEGER(kind=c_int) :: gdalgcpstogeotransform_f

gdalgcpstogeotransform_f = gdalgcpstogeotransform(SIZE(pasgcps), pasgcps(1), padfgeotransform, bapproxok)

END FUNCTION gdalgcpstogeotransform_f


! ========================================
! gdalapplygeotransform
! ========================================
! this unfortunately does not work as ELEMENTAL, padfgeotransform
! should be a scalar derived type
SUBROUTINE gdalapplygeotransform_f_0d(padfgeotransform, &
 dfpixel, dfline, pdfgeox, pdfgeoy)
REAL(kind=c_double),INTENT(in) :: padfgeotransform(6)
REAL(kind=c_double),INTENT(in) :: dfpixel
REAL(kind=c_double),INTENT(in) :: dfline
REAL(kind=c_double),INTENT(out) :: pdfgeox
REAL(kind=c_double),INTENT(out) :: pdfgeoy

pdfgeox = padfgeotransform(1) + &
 dfpixel * padfgeotransform(2) + dfline * padfgeotransform(3)
pdfgeoy = padfgeotransform(4) + &
 dfpixel * padfgeotransform(5) + dfline * padfgeotransform(6)

END SUBROUTINE gdalapplygeotransform_f_0d

SUBROUTINE gdalapplygeotransform_f_1d(padfgeotransform, &
 dfpixel, dfline, pdfgeox, pdfgeoy)
REAL(kind=c_double),INTENT(in) :: padfgeotransform(6)
REAL(kind=c_double),INTENT(in) :: dfpixel(:)
REAL(kind=c_double),INTENT(in) :: dfline(:)
REAL(kind=c_double),INTENT(out) :: pdfgeox(:)
REAL(kind=c_double),INTENT(out) :: pdfgeoy(:)

pdfgeox = padfgeotransform(1) + &
 dfpixel * padfgeotransform(2) + dfline * padfgeotransform(3)
pdfgeoy = padfgeotransform(4) + &
 dfpixel * padfgeotransform(5) + dfline * padfgeotransform(6)

END SUBROUTINE gdalapplygeotransform_f_1d

SUBROUTINE gdalapplygeotransform_f_2d(padfgeotransform, &
 dfpixel, dfline, pdfgeox, pdfgeoy)
REAL(kind=c_double),INTENT(in) :: padfgeotransform(6)
REAL(kind=c_double),INTENT(in) :: dfpixel(:,:)
REAL(kind=c_double),INTENT(in) :: dfline(:,:)
REAL(kind=c_double),INTENT(out) :: pdfgeox(:,:)
REAL(kind=c_double),INTENT(out) :: pdfgeoy(:,:)

pdfgeox = padfgeotransform(1) + &
 dfpixel * padfgeotransform(2) + dfline * padfgeotransform(3)
pdfgeoy = padfgeotransform(4) + &
 dfpixel * padfgeotransform(5) + dfline * padfgeotransform(6)

END SUBROUTINE gdalapplygeotransform_f_2d

SUBROUTINE gdalapplygeotransform_f_3d(padfgeotransform, &
 dfpixel, dfline, pdfgeox, pdfgeoy)
REAL(kind=c_double),INTENT(in) :: padfgeotransform(6)
REAL(kind=c_double),INTENT(in) :: dfpixel(:,:,:)
REAL(kind=c_double),INTENT(in) :: dfline(:,:,:)
REAL(kind=c_double),INTENT(out) :: pdfgeox(:,:,:)
REAL(kind=c_double),INTENT(out) :: pdfgeoy(:,:,:)

pdfgeox = padfgeotransform(1) + &
 dfpixel * padfgeotransform(2) + dfline * padfgeotransform(3)
pdfgeoy = padfgeotransform(4) + &
 dfpixel * padfgeotransform(5) + dfline * padfgeotransform(6)

END SUBROUTINE gdalapplygeotransform_f_3d


! ========================================
! gdaldatasetrasterio
! ========================================
FUNCTION gdaldatasetrasterio_int8(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int8_t),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_int8

FUNCTION gdaldatasetrasterio_int8_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
INTEGER(kind=c_int8_t),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_byte, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_int8_loc


FUNCTION gdaldatasetrasterio_int16(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int16_t),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_int16

FUNCTION gdaldatasetrasterio_int16_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
INTEGER(kind=c_int16_t),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_int16, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_int16_loc


FUNCTION gdaldatasetrasterio_int32(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int32_t),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_int32

FUNCTION gdaldatasetrasterio_int32_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
INTEGER(kind=c_int32_t),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_int32, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_int32_loc


FUNCTION gdaldatasetrasterio_float(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
REAL(kind=c_float),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_float

FUNCTION gdaldatasetrasterio_float_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
REAL(kind=c_float),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_float32, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_float_loc


FUNCTION gdaldatasetrasterio_double(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
REAL(kind=c_double),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_double

FUNCTION gdaldatasetrasterio_double_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
REAL(kind=c_double),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_float64, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_double_loc


FUNCTION gdaldatasetrasterio_float_cmplx(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
COMPLEX(kind=c_float_complex),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_float_cmplx

FUNCTION gdaldatasetrasterio_float_cmplx_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
COMPLEX(kind=c_float_complex),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_cfloat32, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_float_cmplx_loc


FUNCTION gdaldatasetrasterio_double_cmplx(hds, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
COMPLEX(kind=c_double_complex),INTENT(inout) :: pbuffer(:,:,:)
INTEGER(kind=c_int) :: err ! CPLErr

INTEGER(kind=c_int) :: i

err = gdaldatasetrasterio_loc(hds, erwflag, ndsxoff, ndsyoff, &
 SIZE(pbuffer,1), SIZE(pbuffer,2), SIZE(pbuffer,3), pbuffer, &
 (/(i,i=1,SIZE(pbuffer,3))/))

END FUNCTION gdaldatasetrasterio_double_cmplx

FUNCTION gdaldatasetrasterio_double_cmplx_loc(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, nbandcount, pbuffer, panbandcount) result(err)
TYPE(gdaldataseth),VALUE :: hds
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int),INTENT(in) :: nbandcount
COMPLEX(kind=c_double_complex),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize,nbandcount)
INTEGER(kind=c_int),INTENT(in) :: panbandcount(*)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdaldatasetrasterio(hds, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1,1)), &
 ndsxsize, ndsysize, gdt_cfloat64, nbandcount, panbandcount, 0, 0, 0)

END FUNCTION gdaldatasetrasterio_double_cmplx_loc


! ========================================
! gdaldatasetrasterio
! ========================================
FUNCTION gdalrasterio_int8(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int8_t),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_int8

FUNCTION gdalrasterio_int8_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int8_t),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_byte, 0, 0)

END FUNCTION gdalrasterio_int8_loc


FUNCTION gdalrasterio_int16(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int16_t),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_int16

FUNCTION gdalrasterio_int16_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int16_t),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_int16, 0, 0)

END FUNCTION gdalrasterio_int16_loc


FUNCTION gdalrasterio_int32(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int32_t),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_int32

FUNCTION gdalrasterio_int32_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
INTEGER(kind=c_int32_t),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_int32, 0, 0)

END FUNCTION gdalrasterio_int32_loc


FUNCTION gdalrasterio_float(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
REAL(kind=c_float),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_float

FUNCTION gdalrasterio_float_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
REAL(kind=c_float),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_float32, 0, 0)

END FUNCTION gdalrasterio_float_loc


FUNCTION gdalrasterio_double(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
REAL(kind=c_double),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_double

FUNCTION gdalrasterio_double_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
REAL(kind=c_double),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_float64, 0, 0)

END FUNCTION gdalrasterio_double_loc


FUNCTION gdalrasterio_float_cmplx(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
COMPLEX(kind=c_float_complex),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_float_cmplx

FUNCTION gdalrasterio_float_cmplx_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
COMPLEX(kind=c_float_complex),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_cfloat32, 0, 0)

END FUNCTION gdalrasterio_float_cmplx_loc


FUNCTION gdalrasterio_double_cmplx(hband, erwflag, ndsxoff, ndsyoff, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
COMPLEX(kind=c_double_complex),INTENT(inout) :: pbuffer(:,:)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio_loc(hband, erwflag, ndsxoff, ndsyoff, &
   SIZE(pbuffer,1), SIZE(pbuffer,2), pbuffer)

END FUNCTION gdalrasterio_double_cmplx

FUNCTION gdalrasterio_double_cmplx_loc(hband, erwflag, ndsxoff, ndsyoff, ndsxsize, ndsysize, pbuffer) RESULT(err)
TYPE(gdalrasterbandh),VALUE :: hband
INTEGER(kind=c_int),INTENT(in) :: erwflag
INTEGER(kind=c_int),INTENT(in) :: ndsxoff, ndsyoff
INTEGER(kind=c_int),INTENT(in) :: ndsxsize
INTEGER(kind=c_int),INTENT(in) :: ndsysize
COMPLEX(kind=c_double_complex),TARGET,INTENT(inout) :: pbuffer(ndsxsize,ndsysize)
INTEGER(kind=c_int) :: err ! CPLErr

err = gdalrasterio(hband, erwflag, ndsxoff, ndsyoff, &
 ndsxsize, ndsysize, c_loc(pbuffer(1,1)), &
 ndsxsize, ndsysize, gdt_cfloat64, 0, 0)

END FUNCTION gdalrasterio_double_cmplx_loc

SUBROUTINE gdaldatasetbbsize_f(hds, bbxmin, bbymin, bbxmax, bbymax, &
 nx, ny, offsetx, offsety, xmin, ymin, xmax, ymax)
TYPE(gdaldataseth),VALUE :: hds
REAL(kind=c_double),INTENT(in) :: bbxmin
REAL(kind=c_double),INTENT(in) :: bbymin
REAL(kind=c_double),INTENT(in) :: bbxmax
REAL(kind=c_double),INTENT(in) :: bbymax
INTEGER,intent(out) :: nx
INTEGER,intent(out) :: ny
INTEGER,intent(out) :: offsetx
INTEGER,intent(out) :: offsety
REAL(kind=c_double),INTENT(out) :: xmin
REAL(kind=c_double),INTENT(out) :: ymin
REAL(kind=c_double),INTENT(out) :: xmax
REAL(kind=c_double),INTENT(out) :: ymax
!REAL(kind=c_double),INTENT(out) :: sx !< grid step in the x direction
!REAL(kind=c_double),INTENT(out) :: sy !< grid step in the y direction

INTEGER(kind=c_int) :: ier
REAL(kind=c_double) :: geotrans(6), invgeotrans(6), i1r, j1r, i2r, j2r, &
 x1, y1, x2, y2
REAL(kind=c_double),PARAMETER :: epsy = 0.1
INTEGER(kind=c_int) :: i1, j1, i2, j2 !, offsetx, offsety, nx, ny

! ensure (anti)diagonality
ier = gdalgetgeotransform(hds, geotrans)
IF (.NOT.(geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) .AND. &
 .NOT.(geotrans(2) == 0.0_c_double .AND. geotrans(6) == 0.0_c_double)) THEN
  nx = -1
  ny = -1
  RETURN
ENDIF

! compute real indices of bounding box requested
ier = gdalinvgeotransform(geotrans, invgeotrans)
CALL gdalapplygeotransform(invgeotrans, bbxmin, bbymin, i1r, j1r)
CALL gdalapplygeotransform(invgeotrans, bbxmax, bbymax, i2r, j2r)

! compute integer indices of bounding box requested within the domain
i1 = max(nint(min(i1r, i2r) - epsy), 0)
j1 = max(nint(min(j1r, j2r) - epsy), 0)
i2 = min(nint(max(i1r, i2r) + epsy), gdalgetrasterxsize(hds))
j2 = min(nint(max(j1r, j2r) + epsy), gdalgetrasterysize(hds))
offsetx = i1
offsety = j1
nx = max(i2 - i1, 0) ! 0=bounding box outside dataset
ny = max(j2 - j1, 0) ! 0=bounding box outside dataset

! compute output grid corners and steps
CALL gdalapplygeotransform(geotrans, i1 + 0.5_c_double, j1 + 0.5_c_double, &
 x1, y1)
CALL gdalapplygeotransform(geotrans, i2 - 0.5_c_double, j2 - 0.5_c_double, &
 x2, y2)

xmin = min(x1, x2)
ymin = min(y1, y2)
xmax = max(x1, x2)
ymax = max(y1, y2)
!sx = ABS(geotrans(2)) ! improve
!sy = ABS(geotrans(6)) ! improve

END SUBROUTINE gdaldatasetbbsize_f


SUBROUTINE gdaldatasetsimpleread_f(hds, bbxmin, bbymin, bbxmax, bbymax, pbuffer, &
 xmin, ymin, xmax, ymax)
TYPE(gdaldataseth),VALUE :: hds
REAL(kind=c_double),INTENT(in) :: bbxmin
REAL(kind=c_double),INTENT(in) :: bbymin
REAL(kind=c_double),INTENT(in) :: bbxmax
REAL(kind=c_double),INTENT(in) :: bbymax
REAL(kind=c_float),ALLOCATABLE,INTENT(out) :: pbuffer(:,:,:)
REAL(kind=c_double),INTENT(out) :: xmin
REAL(kind=c_double),INTENT(out) :: ymin
REAL(kind=c_double),INTENT(out) :: xmax
REAL(kind=c_double),INTENT(out) :: ymax
!REAL(kind=c_double),INTENT(out) :: sx !< grid step in the x direction
!REAL(kind=c_double),INTENT(out) :: sy !< grid step in the y direction

INTEGER(kind=c_int) :: ier
INTEGER(kind=c_int) :: nx, ny, offsetx, offsety


CALL gdaldatasetbbsize_f(hds, bbxmin, bbymin, bbxmax, bbymax, &
 nx, ny, offsetx, offsety, xmin, ymin, xmax, ymax)
IF (nx < 0 .OR. ny < 0) RETURN ! dataset read error

ALLOCATE(pbuffer(nx, ny, gdalgetrastercount(hds)))
IF (nx == 0 .OR. ny == 0) RETURN ! bounding box outside dataset

ier = gdaldatasetrasterio_f(hds, gf_read, offsetx, offsety, pbuffer)
IF (ier /= 0) THEN
  DEALLOCATE(pbuffer)
  RETURN
ENDIF

! here we should swap/transpose as requested

END SUBROUTINE gdaldatasetsimpleread_f


SUBROUTINE gdalrastersimpleread_f(hband, bbxmin, bbymin, bbxmax, bbymax, pbuffer, &
 xmin, ymin, xmax, ymax)
TYPE(gdalrasterbandh),VALUE :: hband
REAL(kind=c_double),INTENT(in) :: bbxmin
REAL(kind=c_double),INTENT(in) :: bbymin
REAL(kind=c_double),INTENT(in) :: bbxmax
REAL(kind=c_double),INTENT(in) :: bbymax
REAL(kind=c_float),ALLOCATABLE,INTENT(out) :: pbuffer(:,:)
REAL(kind=c_double),INTENT(out) :: xmin
REAL(kind=c_double),INTENT(out) :: ymin
REAL(kind=c_double),INTENT(out) :: xmax
REAL(kind=c_double),INTENT(out) :: ymax
!REAL(kind=c_double),INTENT(out) :: sx !< grid step in the x direction
!REAL(kind=c_double),INTENT(out) :: sy !< grid step in the y direction

INTEGER(kind=c_int) :: ier
INTEGER(kind=c_int) :: nx, ny, offsetx, offsety


CALL gdaldatasetbbsize_f(gdalgetbanddataset(hband), bbxmin, bbymin, bbxmax, bbymax, &
 nx, ny, offsetx, offsety, xmin, ymin, xmax, ymax)
IF (nx < 0 .OR. ny < 0) RETURN ! dataset read error

ALLOCATE(pbuffer(nx, ny))
IF (nx == 0 .OR. ny == 0) RETURN ! bounding box outside dataset

ier = gdalrasterio_f(hband, gf_read, offsetx, offsety, pbuffer)
IF (ier /= 0) THEN
  DEALLOCATE(pbuffer)
  RETURN
ENDIF

! here we should swap/transpose as requested

END SUBROUTINE gdalrastersimpleread_f


FUNCTION gdalmajorobject_fromdataset_new(gdalobject) RESULT(majorobject)
TYPE(gdaldataseth),VALUE :: gdalobject
TYPE(gdalmajorobjecth) :: majorobject
majorobject%ptr = gdalobject%ptr
END FUNCTION gdalmajorobject_fromdataset_new

FUNCTION gdalmajorobject_fromrasterband_new(gdalobject) RESULT(majorobject)
TYPE(gdalrasterbandh),VALUE :: gdalobject
TYPE(gdalmajorobjecth) :: majorobject
majorobject%ptr = gdalobject%ptr
END FUNCTION gdalmajorobject_fromrasterband_new

FUNCTION gdalmajorobject_fromdriver_new(gdalobject) RESULT(majorobject)
TYPE(gdaldriverh),VALUE :: gdalobject
TYPE(gdalmajorobjecth) :: majorobject
majorobject%ptr = gdalobject%ptr
END FUNCTION gdalmajorobject_fromdriver_new

END MODULE gdal