FortranGIS Version 3.0
gdal_test.F90
1PROGRAM gdal_test
2use,INTRINSIC :: iso_c_binding
3USE fortranc
4USE gdal
5IMPLICIT none
6
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(:,:,:)
14
15
16! register all available gdal drivers
17CALL gdalallregister()
18! file name to work on
19file = 'gdal_test.tiff'
20
21! ==== How to create a gdal dataset and export it to a file ====
22
23! get a specific driver (see e.g. gdalinfo --formats)
24print*,'Getting GeoTIFF driver'
25driver = gdalgetdriverbyname('GTiff'//char(0))
26IF (.NOT.gdalassociated(driver)) THEN
27 print*,'Error getting GeoTIFF driver from gdal'
28 stop 1
29ENDIF
30
31! create the dataset with i1xj1 points and 1 raster band of byte data
32print*,'Creating a GeoTIFF gdal dataset'
33i1 = 120
34j1 = 80
35k1 = 3
36! pass a couple of options (BIGTIFF and COMPRESS) to the create function,
37! the strings must have the same legth (possibly compiler-dependent),
38! they will be trimmed by the c_ptr_ptr_new function
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'/))))
41! (/('',i=1,0)/) is a trick to define a zero-length array on the fly,
42! if we do not want to pass any specific option
43!ds = gdalcreate(driver, TRIM(file)//CHAR(0), i1, j1, k1, gdt_byte, &
44! c_ptr_ptr_getobject(c_ptr_ptr_new((/('',i=1,0)/))))
45IF (.NOT.gdalassociated(ds)) THEN
46 print*,'Error creating a GeoTIFF dataset on file ',trim(file)
47 stop 1
48ENDIF
49
50! this seems to be useless here, depends on file type
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)
55
56! create a dummy array of real data, we stick to integer <= 255 in
57! order to fit in a byte and not to lose precision in read test
58ALLOCATE(z3(i1,j1,k1))
59DO k = 1, k1
60 DO j = 1, j1
61 DO i = 1, i1
62 z3(i,j,k) = i/2+j*mod(k,2)+(j1-j)*(1-mod(k,2))
63 ENDDO
64 ENDDO
65ENDDO
66
67! write data to dataset with the simplified Fortran interface, one
68! raster band is written starting from the upper(?) left corner, the
69! conversion from the type of z3 (real) to the gdt_byte type specified
70! with gdalcreate is done internally by gdal
71print*,'Writing data to dataset'
72print*,'using the simplified Fortran interface'
73ierr = gdaldatasetrasterio_f(ds, gf_write, 0, 0, z3)
74IF (ierr /= 0) THEN
75 print*,'Error writing data to GeoTIFF dataset on file ',trim(file)
76 stop 1
77ENDIF
78
79CALL gdalclose(ds)
80
81! ==== How to read a gdal dataset from a file, simplified Fortran interface ====
82
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)
88 stop 1
89ENDIF
90
91print*,'Getting the affine transformation'
92ierr = gdalgetgeotransform(ds, gt)
93! an error is acceptable since no transformation had been defined
94!IF (ierr /= 0) THEN
95! PRINT*,'Error getting the affine transformation from dataset on file ',TRIM(file)
96! STOP 1
97!ENDIF
98print*,'The affine transformation matrix is: ',gt
99
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
106 stop 1
107ENDIF
108print*,'The x/y size of the raster is: ',i2,j2
109
110IF (k1 /= k2) THEN
111 print*,'Error, wrong raster band number ',k1,k2
112 stop 1
113ENDIF
114print*,'The number of raster bands is: ',k2
115
116! apply the affine transformation to some coordinates
117! original gdal version
118CALL gdalapplygeotransform(gt, 0.5_c_double, 0.5_c_double, x1, y1)
119! Fortran elemental version
120CALL gdalapplygeotransform_f(gt, 0.5_c_double, 0.5_c_double, x2, y2)
121IF (x1 /= x2 .OR. y1 /= y2) THEN ! this check should be relaxed
122 print*,'Error gdal and Fortran version of GDALApplyGeoTransform'
123 print*,'give different results'
124 stop 1
125ENDIF
126
127CALL gdalapplygeotransform_f(gt, 0.5_c_double, 0.5_c_double, x1, y1)
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
132
133ALLOCATE(z(i2,j2))
134
135! get the first raster band
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)
140 stop 1
141ENDIF
142
143! read data from first raster band with the simplified Fortran interface,
144! raster band is read starting from the upper(?) left corner
145print*,'Reading data from dataset'
146ierr = gdalrasterio_f(band, gf_read, 0, 0, z)
147IF (ierr /= 0) THEN
148 print*,'Error reading data from GeoTIFF dataset on file ',trim(file)
149 print*,'with simplified Fortran interface'
150 stop 1
151ENDIF
152
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)
156
157CALL gdalclose(ds)
158
159! ==== How to read a gdal dataset from a file, even more simplified Fortran interface ====
160
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)
166 stop 1
167ENDIF
168
169! read data from first raster band with the even more simplified Fortran interface,
170! raster band is read starting from the upper(?) left corner
171print*,'Reading data from dataset'
172CALL gdaldatasetsimpleread_f(ds, 5._c_double, 5._c_double, 20._c_double, 15._c_double, &
173 zr3, x1, y1, x2, y2)
174
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'
178 stop 1
179ENDIF
180
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
184
185CALL gdalclose(ds)
186
187DEALLOCATE(z3,z,zr3)
188
189END PROGRAM gdal_test
190
Constructor for a c_ptr_ptr object.
Definition fortranc.F90:185
Interface to a Fortran version of gdalapplygeotransform working on scalars, 1-d, 2-d and 3-d arrays.
Definition gdal.F90:314
Simplified Fortran generic interface to the gdaldatasetrasterio C function.
Definition gdal.F90:343
Simplified Fortran generic interface to the gdalrasterio C function.
Definition gdal.F90:372
Utility module for supporting Fortran 2003 C language interface module.
Definition fortranc.F90:103
Fortran 2003 interface to the gdal http://www.gdal.org/ library.
Definition gdal.F90:201