gdal_test.F90

PROGRAM gdal_test
use,INTRINSIC :: iso_c_binding
USE fortranc
USE gdal
IMPLICIT none

TYPE(gdaldriverh) :: driver
TYPE(gdaldataseth) :: ds
TYPE(gdalrasterbandh) :: band
CHARACTER(len=512) :: file
REAL(kind=c_double) :: x1, y1, x2, y2, gt(6)
INTEGER(kind=c_int) :: i1, j1, k1, i2, j2, k2, i, j, k, ierr
REAL,ALLOCATABLE :: z(:,:), z3(:,:,:), zr3(:,:,:)


! register all available gdal drivers
CALL gdalallregister()
! file name to work on
file = 'gdal_test.tiff'

! ==== How to create a gdal dataset and export it to a file ====

! get a specific driver (see e.g. gdalinfo --formats)
print*,'Getting GeoTIFF driver'
driver = gdalgetdriverbyname('GTiff'//char(0))
IF (.NOT.gdalassociated(driver)) THEN
  print*,'Error getting GeoTIFF driver from gdal'
  stop 1
ENDIF

! create the dataset with i1xj1 points and 1 raster band of byte data
print*,'Creating a GeoTIFF gdal dataset'
i1 = 120
j1 = 80
k1 = 3
ds = gdalcreate(driver, trim(file)//char(0), i1, j1, k1, gdt_byte, &
 c_ptr_ptr_getobject(c_ptr_ptr_new((/('',i=1,0)/))))
! (/('',i=1,0)/) is a trick to define a zero-length array on the fly,
! since we do not want to pass any specific option
IF (.NOT.gdalassociated(ds)) THEN
  print*,'Error creating a GeoTIFF dataset on file ',trim(file)
  stop 1
ENDIF

! this seems to be useless here, depends on file type
print*,'Setting color interpretation to RGB'
ierr =  gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 1), gci_redband)
ierr =  gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 2), gci_greenband)
ierr =  gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 3), gci_blueband)

! create a dummy array of real data, we stick to integer <= 255 in
! order to fit in a byte and not to lose precision in read test
ALLOCATE(z3(i1,j1,k1))
DO k = 1, k1
  DO j = 1, j1
    DO i = 1, i1
      z3(i,j,k) = i/2+j*mod(k,2)+(j1-j)*(1-mod(k,2))
    ENDDO
  ENDDO
ENDDO

! write data to dataset with the simplified Fortran interface, one
! raster band is written starting from the upper(?) left corner, the
! conversion from the type of z3 (real) to the gdt_byte type specified
! with gdalcreate is done internally by gdal
print*,'Writing data to dataset'
print*,'using the simplified Fortran interface'
ierr = gdaldatasetrasterio_f(ds, gf_write, 0, 0, z3)
IF (ierr /= 0) THEN
  print*,'Error writing data to GeoTIFF dataset on file ',trim(file)
  stop 1
ENDIF

CALL gdalclose(ds)

! ==== How to read a gdal dataset from a file, simplified Fortran interface ====

print*,'Opening a GeoTIFF gdal dataset for reading'
print*,'using the simplified Fortran interface'
ds = gdalopen(trim(file)//char(0), ga_readonly)
IF (.NOT.gdalassociated(ds)) THEN
  print*,'Error opening dataset on file ',trim(file)
  stop 1
ENDIF

print*,'Getting the affine transformation'
ierr = gdalgetgeotransform(ds, gt)
! an error is acceptable since no transformation had been defined
!IF (ierr /= 0) THEN
!  PRINT*,'Error getting the affine transformation from dataset on file ',TRIM(file)
!  STOP 1
!ENDIF
print*,'The affine transformation matrix is: ',gt

print*,'Getting the dataset size'
i2 = gdalgetrasterxsize(ds)
j2 = gdalgetrasterysize(ds)
k2 = gdalgetrastercount(ds)
IF (i1 /= i2 .OR. j1 /= j2) THEN
  print*,'Error, wrong dataset x or y size ',i1,i2,j1,j2
  stop 1
ENDIF
print*,'The x/y size of the raster is: ',i2,j2

IF (k1 /= k2) THEN
  print*,'Error, wrong raster band number ',k1,k2
  stop 1
ENDIF
print*,'The number of raster bands is: ',k2

! apply the affine transformation to some coordinates
! original gdal version
CALL gdalapplygeotransform(gt, 0.5_c_double, 0.5_c_double, x1, y1)
! Fortran elemental version
CALL gdalapplygeotransform_f(gt, 0.5_c_double, 0.5_c_double, x2, y2)
IF (x1 /= x2 .OR. y1 /= y2) THEN ! this check should be relaxed
  print*,'Error gdal and Fortran version of GDALApplyGeoTransform'
  print*,'give different results'
  stop 1
ENDIF

CALL gdalapplygeotransform_f(gt, 0.5_c_double, 0.5_c_double, x1, y1)
CALL gdalapplygeotransform_f(gt, &
 REAL(i1, kind=c_double)-0.5_c_double, &
 REAL(j1, kind=c_double)-0.5_c_double, x2, y2)
print*,'The raster coordinates of the corners are: ',x1,y1,x2,y2

ALLOCATE(z(i2,j2))

! get the first raster band
print*,'Getting a raster band'
band = gdalgetrasterband(ds, 1)
IF (.NOT.gdalassociated(band)) THEN
  print*,'Error getting raster band from GeoTIFF dataset on file ',trim(file)
  stop 1
ENDIF

! read data from first raster band with the simplified Fortran interface,
! raster band is read starting from the upper(?) left corner
print*,'Reading data from dataset'
ierr = gdalrasterio_f(band, gf_read, 0, 0, z)
IF (ierr /= 0) THEN
  print*,'Error reading data from GeoTIFF dataset on file ',trim(file)
  print*,'with simplified Fortran interface'
  stop 1
ENDIF

print*,'The sum of the buffer read is: ',sum(z)
print*,'Average error after write/read process: ',&
 sum(abs(z(:,:) - z3(:,:,1)))/(i2*j2)

CALL gdalclose(ds)

! ==== How to read a gdal dataset from a file, even more simplified Fortran interface ====

print*,'Opening a GeoTIFF gdal dataset for reading'
print*,'using the even more simplified Fortran interface'
ds = gdalopen(trim(file)//char(0), ga_readonly)
IF (.NOT.gdalassociated(ds)) THEN
  print*,'Error opening dataset on file ',trim(file)
  stop 1
ENDIF

! read data from first raster band with the even more simplified Fortran interface,
! raster band is read starting from the upper(?) left corner
print*,'Reading data from dataset'
CALL gdaldatasetsimpleread_f(ds, 5._c_double, 5._c_double, 20._c_double, 15._c_double, &
 zr3, x1, y1, x2, y2)

IF (.NOT.ALLOCATED(zr3)) THEN
  print*,'Error reading data from GeoTIFF dataset on file ',trim(file)
  print*,'with even more simplified Fortran interface'
  stop 1
ENDIF

print*,'The shape of the buffer read is: ',shape(zr3)
print*,'The sum of the buffer read is: ',sum(zr3)
print*,'The raster coordinates of the corners are: ',x1,y1,x2,y2

CALL gdalclose(ds)

DEALLOCATE(z3,z,zr3)

END PROGRAM gdal_test