FortranGIS  Version 3.0
gdal_test.F90
1 PROGRAM gdal_test
2 use,INTRINSIC :: iso_c_binding
3 USE fortranc
4 USE gdal
5 IMPLICIT none
6 
7 TYPE(gdaldriverh) :: driver
8 TYPE(gdaldataseth) :: ds
9 TYPE(gdalrasterbandh) :: band
10 CHARACTER(len=512) :: file
11 REAL(kind=c_double) :: x1, y1, x2, y2, gt(6)
12 INTEGER(kind=c_int) :: i1, j1, k1, i2, j2, k2, i, j, k, ierr
13 REAL,ALLOCATABLE :: z(:,:), z3(:,:,:), zr3(:,:,:)
14 
15 
16 ! register all available gdal drivers
17 CALL gdalallregister()
18 ! file name to work on
19 file = '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)
24 print*,'Getting GeoTIFF driver'
25 driver = gdalgetdriverbyname('GTiff'//char(0))
26 IF (.NOT.gdalassociated(driver)) THEN
27  print*,'Error getting GeoTIFF driver from gdal'
28  stop 1
29 ENDIF
30 
31 ! create the dataset with i1xj1 points and 1 raster band of byte data
32 print*,'Creating a GeoTIFF gdal dataset'
33 i1 = 120
34 j1 = 80
35 k1 = 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
39 ds = 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)/))))
45 IF (.NOT.gdalassociated(ds)) THEN
46  print*,'Error creating a GeoTIFF dataset on file ',trim(file)
47  stop 1
48 ENDIF
49 
50 ! this seems to be useless here, depends on file type
51 print*,'Setting color interpretation to RGB'
52 ierr = gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 1), gci_redband)
53 ierr = gdalsetrastercolorinterpretation(gdalgetrasterband(ds, 2), gci_greenband)
54 ierr = 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
58 ALLOCATE(z3(i1,j1,k1))
59 DO 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
65 ENDDO
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
71 print*,'Writing data to dataset'
72 print*,'using the simplified Fortran interface'
73 ierr = gdaldatasetrasterio_f(ds, gf_write, 0, 0, z3)
74 IF (ierr /= 0) THEN
75  print*,'Error writing data to GeoTIFF dataset on file ',trim(file)
76  stop 1
77 ENDIF
78 
79 CALL gdalclose(ds)
80 
81 ! ==== How to read a gdal dataset from a file, simplified Fortran interface ====
82 
83 print*,'Opening a GeoTIFF gdal dataset for reading'
84 print*,'using the simplified Fortran interface'
85 ds = gdalopen(trim(file)//char(0), ga_readonly)
86 IF (.NOT.gdalassociated(ds)) THEN
87  print*,'Error opening dataset on file ',trim(file)
88  stop 1
89 ENDIF
90 
91 print*,'Getting the affine transformation'
92 ierr = 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
98 print*,'The affine transformation matrix is: ',gt
99 
100 print*,'Getting the dataset size'
101 i2 = gdalgetrasterxsize(ds)
102 j2 = gdalgetrasterysize(ds)
103 k2 = gdalgetrastercount(ds)
104 IF (i1 /= i2 .OR. j1 /= j2) THEN
105  print*,'Error, wrong dataset x or y size ',i1,i2,j1,j2
106  stop 1
107 ENDIF
108 print*,'The x/y size of the raster is: ',i2,j2
109 
110 IF (k1 /= k2) THEN
111  print*,'Error, wrong raster band number ',k1,k2
112  stop 1
113 ENDIF
114 print*,'The number of raster bands is: ',k2
115 
116 ! apply the affine transformation to some coordinates
117 ! original gdal version
118 CALL gdalapplygeotransform(gt, 0.5_c_double, 0.5_c_double, x1, y1)
119 ! Fortran elemental version
120 CALL gdalapplygeotransform_f(gt, 0.5_c_double, 0.5_c_double, x2, y2)
121 IF (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
125 ENDIF
126 
127 CALL gdalapplygeotransform_f(gt, 0.5_c_double, 0.5_c_double, x1, y1)
128 CALL gdalapplygeotransform_f(gt, &
129  REAL(i1, kind=c_double)-0.5_c_double, &
130  REAL(j1, kind=c_double)-0.5_c_double, x2, y2)
131 print*,'The raster coordinates of the corners are: ',x1,y1,x2,y2
132 
133 ALLOCATE(z(i2,j2))
134 
135 ! get the first raster band
136 print*,'Getting a raster band'
137 band = gdalgetrasterband(ds, 1)
138 IF (.NOT.gdalassociated(band)) THEN
139  print*,'Error getting raster band from GeoTIFF dataset on file ',trim(file)
140  stop 1
141 ENDIF
142 
143 ! read data from first raster band with the simplified Fortran interface,
144 ! raster band is read starting from the upper(?) left corner
145 print*,'Reading data from dataset'
146 ierr = gdalrasterio_f(band, gf_read, 0, 0, z)
147 IF (ierr /= 0) THEN
148  print*,'Error reading data from GeoTIFF dataset on file ',trim(file)
149  print*,'with simplified Fortran interface'
150  stop 1
151 ENDIF
152 
153 print*,'The sum of the buffer read is: ',sum(z)
154 print*,'Average error after write/read process: ',&
155  sum(abs(z(:,:) - z3(:,:,1)))/(i2*j2)
156 
157 CALL gdalclose(ds)
158 
159 ! ==== How to read a gdal dataset from a file, even more simplified Fortran interface ====
160 
161 print*,'Opening a GeoTIFF gdal dataset for reading'
162 print*,'using the even more simplified Fortran interface'
163 ds = gdalopen(trim(file)//char(0), ga_readonly)
164 IF (.NOT.gdalassociated(ds)) THEN
165  print*,'Error opening dataset on file ',trim(file)
166  stop 1
167 ENDIF
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
171 print*,'Reading data from dataset'
172 CALL gdaldatasetsimpleread_f(ds, 5._c_double, 5._c_double, 20._c_double, 15._c_double, &
173  zr3, x1, y1, x2, y2)
174 
175 IF (.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
179 ENDIF
180 
181 print*,'The shape of the buffer read is: ',shape(zr3)
182 print*,'The sum of the buffer read is: ',sum(zr3)
183 print*,'The raster coordinates of the corners are: ',x1,y1,x2,y2
184 
185 CALL gdalclose(ds)
186 
187 DEALLOCATE(z3,z,zr3)
188 
189 END 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