2use,
INTRINSIC :: iso_c_binding
7TYPE(gdaldriverh) :: driver
8TYPE(gdaldataseth) :: ds
9TYPE(gdalrasterbandh) :: band
10CHARACTER(len=512) :: file
11REAL(kind=c_double) :: x1, y1, x2, y2, gt(6)
12INTEGER(kind=c_int) :: i1, j1, k1, i2, j2, k2, i, j, k, ierr
13REAL,
ALLOCATABLE :: z(:,:), z3(:,:,:), zr3(:,:,:)
19file =
'gdal_test.tiff'
24print*,
'Getting GeoTIFF driver'
25driver = gdalgetdriverbyname(
'GTiff'//char(0))
26IF (.NOT.gdalassociated(driver))
THEN
27 print*,
'Error getting GeoTIFF driver from gdal'
32print*,
'Creating a GeoTIFF gdal dataset'
39ds = gdalcreate(driver, trim(file)//char(0), i1, j1, k1, gdt_byte, &
40 c_ptr_ptr_getobject(
c_ptr_ptr_new((/
'BIGTIFF=YES ',
'COMPRESS=DEFLATE'/))))
45IF (.NOT.gdalassociated(ds))
THEN
46 print*,
'Error creating a GeoTIFF dataset on file ',trim(file)
51print*,
'Setting color interpretation to RGB'
52ierr = gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 1), gci_redband)
53ierr = gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 2), gci_greenband)
54ierr = gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 3), gci_blueband)
62 z3(i,j,k) = i/2+j*mod(k,2)+(j1-j)*(1-mod(k,2))
71print*,
'Writing data to dataset'
72print*,
'using the simplified Fortran interface'
75 print*,
'Error writing data to GeoTIFF dataset on file ',trim(file)
83print*,
'Opening a GeoTIFF gdal dataset for reading'
84print*,
'using the simplified Fortran interface'
85ds = gdalopen(trim(file)//char(0), ga_readonly)
86IF (.NOT.gdalassociated(ds))
THEN
87 print*,
'Error opening dataset on file ',trim(file)
91print*,
'Getting the affine transformation'
92ierr = gdalgetgeotransform(ds, gt)
98print*,
'The affine transformation matrix is: ',gt
100print*,
'Getting the dataset size'
101i2 = gdalgetrasterxsize(ds)
102j2 = gdalgetrasterysize(ds)
103k2 = gdalgetrastercount(ds)
104IF (i1 /= i2 .OR. j1 /= j2)
THEN
105 print*,
'Error, wrong dataset x or y size ',i1,i2,j1,j2
108print*,
'The x/y size of the raster is: ',i2,j2
111 print*,
'Error, wrong raster band number ',k1,k2
114print*,
'The number of raster bands is: ',k2
118CALL gdalapplygeotransform(gt, 0.5_c_double, 0.5_c_double, x1, y1)
121IF (x1 /= x2 .OR. y1 /= y2)
THEN
122 print*,
'Error gdal and Fortran version of GDALApplyGeoTransform'
123 print*,
'give different results'
129 REAL(i1, kind=c_double)-0.5_c_double, &
130 REAL(j1, kind=c_double)-0.5_c_double, x2, y2)
131print*,
'The raster coordinates of the corners are: ',x1,y1,x2,y2
136print*,
'Getting a raster band'
137band = gdalgetrasterband(ds, 1)
138IF (.NOT.gdalassociated(band))
THEN
139 print*,
'Error getting raster band from GeoTIFF dataset on file ',trim(file)
145print*,
'Reading data from dataset'
148 print*,
'Error reading data from GeoTIFF dataset on file ',trim(file)
149 print*,
'with simplified Fortran interface'
153print*,
'The sum of the buffer read is: ',sum(z)
154print*,
'Average error after write/read process: ',&
155 sum(abs(z(:,:) - z3(:,:,1)))/(i2*j2)
161print*,
'Opening a GeoTIFF gdal dataset for reading'
162print*,
'using the even more simplified Fortran interface'
163ds = gdalopen(trim(file)//char(0), ga_readonly)
164IF (.NOT.gdalassociated(ds))
THEN
165 print*,
'Error opening dataset on file ',trim(file)
171print*,
'Reading data from dataset'
172CALL gdaldatasetsimpleread_f(ds, 5._c_double, 5._c_double, 20._c_double, 15._c_double, &
175IF (.NOT.
ALLOCATED(zr3))
THEN
176 print*,
'Error reading data from GeoTIFF dataset on file ',trim(file)
177 print*,
'with even more simplified Fortran interface'
181print*,
'The shape of the buffer read is: ',shape(zr3)
182print*,
'The sum of the buffer read is: ',sum(zr3)
183print*,
'The raster coordinates of the corners are: ',x1,y1,x2,y2
Constructor for a c_ptr_ptr object.
Simplified Fortran generic interface to the gdaldatasetrasterio C function.
Simplified Fortran generic interface to the gdalrasterio C function.
Utility module for supporting Fortran 2003 C language interface module.
Fortran 2003 interface to the gdal http://www.gdal.org/ library.