FortranGIS Version 3.0
shapelib_test.F90
1! Copyright 2011 Davide Cesari <dcesari69 at gmail dot com>
2!
3! This file is part of FortranGIS.
4!
5! FortranGIS is free software: you can redistribute it and/or modify
6! it under the terms of the GNU Lesser General Public License as
7! published by the Free Software Foundation, either version 3 of the
8! License, or (at your option) any later version.
9!
10! FortranGIS is distributed in the hope that it will be useful, but
11! WITHOUT ANY WARRANTY; without even the implied warranty of
12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13! Lesser General Public License for more details.
14!
15! You should have received a copy of the GNU Lesser General Public
16! License along with FortranGIS. If not, see
17! <http://www.gnu.org/licenses/>.
18
19PROGRAM shapelib_test
20use,INTRINSIC :: iso_c_binding
21USE shapelib
22IMPLICIT NONE
23
24INTEGER,PARAMETER :: lencharattr=40, nshp=4, tshp=shpt_polygonz
25
26TYPE(shpfileobject) :: shphandle
27TYPE(shpobject) :: shpobj
28INTEGER :: i, j
29CHARACTER(len=1024) :: filename
30
31INTEGER :: nshpr, tshpr, nfield, nrec, nd
32REAL(kind=c_double) :: minbound(4), maxbound(4)
33CHARACTER(len=lencharattr) :: charattrr
34INTEGER :: intattrr
35REAL(kind=c_double) :: doubleattrr
36
37!CALL get_command_argument(1,filename)
38!IF (filename == '') THEN
39! PRINT'(A)','Usage: shape_test <shp_file>'
40! STOP
41!ENDIF
42filename = 'testshape'
43
44! ==== How to create a shapefile ====
45
46! create a new shapefile object with data of type tshp (polygon)
47! and associate it to a file, filename does not include extension
48shphandle = shpcreate(trim(filename), tshp)
49! error check
50IF (shpfileisnull(shphandle) .OR. dbffileisnull(shphandle)) THEN
51 print*,'Error opening ',trim(filename),' for writing'
52 stop 1
53ENDIF
54
55! add 3 dbf field, of type character, integer and double respectively
56j = dbfaddfield(shphandle, 'name', ftstring, lencharattr, 0)
57IF (j /= 0) THEN
58 print*,'Error in dbfaddfield',0,j
59 stop 1
60ENDIF
61j = dbfaddfield(shphandle, 'number', ftinteger, 10, 0)
62IF (j /= 1) THEN
63 print*,'Error in dbfaddfield',1,j
64 stop 1
65ENDIF
66j = dbfaddfield(shphandle, 'size', ftdouble, 30, 20)
67IF (j /= 2) THEN
68 print*,'Error in dbfaddfield',2,j
69 stop 1
70ENDIF
71
72! add nshp shapes, of different lengths
73DO i = 0, nshp - 1
74 print*,'Creating shape',i
75! create a shape object with the "simple" method
76! for each shape 3 components are added x, y, z
77! the type of shape has to be repeated here
78! makesimpleshp() is an utility function returning an array
79 shpobj = shpcreatesimpleobject(tshp, &
80 SIZE(makesimpleshp(i, 0)), &
81 makesimpleshp(i, 0), &
82 makesimpleshp(i, 1), &
83 makesimpleshp(i, 2))
84
85! write the shape object to the shapefile object as i-th element
86! -1 = append
87 j = shpwriteobject(shphandle, -1, shpobj)
88 IF (j /= i) THEN
89 print*,'Error in shpwriteobject',i,j
90 stop 1
91 ENDIF
92
93! destroy the shape object to avoid memory leaks
94 CALL shpdestroyobject(shpobj)
95
96! write the 3 dbf attributes of different types for the i-th shape
97! object to the shapefile object
98! makechardbf(), makeintdbf() and makedoubledbf() are utility functions
99! returning an attribute of the proper type
100 j = dbfwriteattribute(shphandle, i, 0, makechardbf(i))
101 IF (j /= 1) THEN
102 print*,'Error in dbfwriteattribute, char',j
103 stop 1
104 ENDIF
105 j = dbfwriteattribute(shphandle, i, 1, makeintdbf(i))
106 IF (j /= 1) THEN
107 print*,'Error in dbfwriteattribute, int',j
108 stop 1
109 ENDIF
110 j = dbfwriteattribute(shphandle, i, 2, makedoubledbf(i))
111 IF (j /= 1) THEN
112 print*,'Warning in dbfwriteattribute, double',j
113 ENDIF
114
115ENDDO
116
117! close the shapefile object
118CALL shpclose(shphandle)
119
120
121! ==== How to read a shapefile ====
122
123! open an exixting shapefile and associate it to a shapefile object
124! filename does not include extension
125shphandle = shpopen(trim(filename), 'rb')
126! error check
127IF (shpfileisnull(shphandle) .OR. dbffileisnull(shphandle)) THEN
128 print*,'Error opening ',trim(filename),' for reading'
129 stop 1
130ENDIF
131
132! get general information about the shapefile object
133CALL shpgetinfo(shphandle, nshpr, tshpr, minbound, maxbound, nfield, nrec)
134IF (nshpr /= nshp) THEN
135 print*,'Error in shpgetinfo, wrong number of shapes',nshp,nshpr
136 stop 1
137ENDIF
138IF (tshpr /= tshp) THEN
139 print*,'Error in shpgetinfo, wrong type of shapes',tshp,tshpr
140 stop 1
141ENDIF
142IF (nfield /= 3) THEN
143 print*,'Error in shpgetinfo, wrong number of fields',3,nfield
144 stop 1
145ENDIF
146IF (nrec /= nshp) THEN
147 print*,'Error in shpgetinfo, wrong number of records',nshp,nrec
148 stop 1
149ENDIF
150
151! read the nshp shapes
152DO i = 0, nshp - 1
153 print*,'Checking shape',i
154! read the i-th shape from the shapefile object and obtain a shape object
155 shpobj = shpreadobject(shphandle, i)
156! error check
157 IF (shpisnull(shpobj)) THEN
158 print*,'Error in shpreadobject',i
159 stop 1
160 ENDIF
161
162! now access all the components of the shape object
163! number of vertices
164 IF (shpobj%nvertices /= SIZE(makesimpleshp(i,0))) THEN
165 print*,'Error in shpreadobject, wrong number of vertices',i,&
166 SIZE(makesimpleshp(i,0)),shpobj%nvertices
167 stop 1
168 ENDIF
169! x shape data
170 IF (any(shpobj%padfx(:) /= makesimpleshp(i,0))) THEN
171 print*,'Error in shpreadobject, discrepancies in x',i
172 print*,makesimpleshp(i,0)
173 print*,shpobj%padfx(:)
174 stop 1
175 ENDIF
176! y shape data
177 IF (any(shpobj%padfy(:) /= makesimpleshp(i,1))) THEN
178 print*,'Error in shpreadobject, discrepancies in y',i
179 print*,makesimpleshp(i,1)
180 print*,shpobj%padfy(:)
181 stop 1
182 ENDIF
183! z shape data
184 IF (any(shpobj%padfz(:) /= makesimpleshp(i,2))) THEN
185 print*,'Error in shpreadobject, discrepancies in z',i
186 print*,makesimpleshp(i,2)
187 print*,shpobj%padfz(:)
188 stop 1
189 ENDIF
190
191! destroy the shape object to avoid memory leaks
192! notice that for accessing dbf attributes the shape object is not required
193 CALL shpdestroyobject(shpobj)
194
195! get the dbf attributes for the i-th shape in the shapefile object
196! first field (character)
197 CALL dbfreadattribute(shphandle, i, 0, charattrr)
198 IF (charattrr /= makechardbf(i)) THEN
199 print*,'Error in dbfreadattribute, discrepancies in char'
200 print*,makechardbf(i)
201 print*,charattrr
202 stop 1
203 ENDIF
204! second field (integer)
205 CALL dbfreadattribute(shphandle, i, 1, intattrr)
206 IF (intattrr /= makeintdbf(i)) THEN
207 print*,'Error in dbfreadattribute, discrepancies in int'
208 print*,makeintdbf(i)
209 print*,intattrr
210 stop 1
211 ENDIF
212! third field (double)
213 CALL dbfreadattribute(shphandle, i, 2, doubleattrr)
214 IF (doubleattrr /= makedoubledbf(i)) THEN
215! here discreapancies are tolerated
216 print*,'Warning in dbfreadattribute, discrepancies in double'
217 print*,makedoubledbf(i)
218 print*,doubleattrr
219 ENDIF
220
221ENDDO
222
223! test whether an attribute is null
224IF (dbfisattributenull(shphandle, 0, 0)) THEN
225 print*,'Error in dbfisattributenull, non null attribute returned null'
226 stop 1
227ENDIF
228IF (.NOT.dbfisattributenull(shphandle, 3, 0)) THEN
229 print*,'Error in dbfisattributenull, null attribute returned non null'
230 stop 1
231ENDIF
232
233! close the shapefile object
234CALL shpclose(shphandle)
235
236CONTAINS
237
238! Functions for generating predictable shp and dbf values
239FUNCTION makesimpleshp(nshp, ncoord) RESULT(shp)
240INTEGER,INTENT(in) :: nshp, ncoord
241REAL(kind=c_double) :: shp(nshp+2)
242
243INTEGER :: i
244
245shp(:) = (/(-100.0_c_double + &
246 10.0_c_double*i + 100.0_c_double*nshp + 1000.0_c_double*ncoord, &
247 i=1, SIZE(shp))/)
248
249END FUNCTION makesimpleshp
250
251FUNCTION makechardbf(nshp) RESULT(dbf)
252INTEGER,INTENT(in) :: nshp
253CHARACTER(len=lencharattr) :: dbf
254
255INTEGER :: i
256
257IF (nshp == 3) THEN
258 dbf= ' '
259ELSE
260 DO i = 1, len(dbf)
261 dbf(i:i) = char(32 + mod(i+2*nshp,32))
262 ENDDO
263ENDIF
264
265
266END FUNCTION makechardbf
267
268FUNCTION makeintdbf(nshp) RESULT(dbf)
269INTEGER,INTENT(in) :: nshp
270INTEGER :: dbf
271
272dbf = -118 + 47*nshp
273
274END FUNCTION makeintdbf
275
276FUNCTION makedoubledbf(nshp) RESULT(dbf)
277INTEGER,INTENT(in) :: nshp
278REAL(kind=c_double) :: dbf
279
280dbf = -5.894823e+12_c_double + 8.4827943e+11*nshp
281
282END FUNCTION makedoubledbf
283
284END PROGRAM shapelib_test
285
Interface to SUBROUTINEs for reading dbf attributes.
Definition: shapelib.F90:138
Interface to FUNCTIONs for setting dbf attributes.
Definition: shapelib.F90:156
Fortran 2003 interface to the shapelib http://shapelib.maptools.org/ library.
Definition: shapelib.F90:53