FortranGIS Version 3.0
shp_fxtr.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 shp_fxtr
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) :: filein, fileout, filereg
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
37CALL get_command_argument(1,filein)
38CALL get_command_argument(2,fileout)
39CALL get_command_argument(3,filereg)
40IF (filein == '' .OR. fileout == '' .OR. filereg == '') THEN
41 print'(A)','Convert a shapefile (with polygons) in the dump format required by'
42 print'(A)','the shape2fxtr.pl utilities of fieldextra and directly in the'
43 print'(A)','region file format required by fieldextra.'
44 print'(A)','Usage: shp_fxtr <shp_file> <dump_file> <reg_file>'
45 stop
46ENDIF
47
48! open an exixting shapefile and associate it to a shapefile object
49! filename does not include extension
50shphandle = shpopen(trim(filein), 'rb')
51! error check
52IF (shpfileisnull(shphandle) .OR. dbffileisnull(shphandle)) THEN
53 print*,'Error opening ',trim(filein),' for reading'
54 stop 1
55ENDIF
56
57! get general information about the shapefile object
58CALL shpgetinfo(shphandle, nshpr, tshpr, minbound, maxbound, nfield, nrec)
59print*,'number and type of shapes:',nshpr,tshpr
60
61OPEN(10, file=fileout)
62OPEN(11, file=filereg)
63! read the nshp shapes
64DO i = 0, nshpr - 1
65 print*,'Checking shape',i
66! read the i-th shape from the shapefile object and obtain a shape object
67 shpobj = shpreadobject(shphandle, i)
68! error check
69 IF (shpisnull(shpobj)) THEN
70 print*,'Error in shpreadobject',i
71 stop 1
72 ENDIF
73
74! now access all the components of the shape object
75! number of vertices
76 print*,'nvertices:',shpobj%nvertices
77 IF (ASSOCIATED(shpobj%panpartstart)) THEN
78! PRINT*,'nparts:',SIZE(shpobj%panpartstart)
79! PRINT*,shpobj%panpartstart
80 IF (SIZE(shpobj%panpartstart) > 1) THEN
81 print*,'Warning, multipart shapes not supported in the conversion'
82 print*,'skipping shape ',i+1
83 cycle
84 ENDIF
85 ENDIF
86 IF (shpobj%nvertices >= 3) THEN
87! simple shp dump file
88 DO j = 1, shpobj%nvertices
89 WRITE(10,'(F0.5,1X,F0.5,1X,I0)')shpobj%padfx(j),shpobj%padfy(j),i+1
90 ENDDO
91! region file
92 WRITE(11,'(I4,1X,I6,1X,''region_'',I0)')0,shpobj%nvertices,i+1
93 WRITE(11,'(*(E10.5,1X))')shpobj%padfy(:) ! warning, f2008-style comment
94 WRITE(11,'(*(E10.5,1X))')shpobj%padfx(:)
95 WRITE(11,'()')
96 ENDIF
97
98! destroy the shape object to avoid memory leaks
99! notice that for accessing dbf attributes the shape object is not required
100 CALL shpdestroyobject(shpobj)
101
102ENDDO
103
104! close the shapefile object
105CALL shpclose(shphandle)
106CLOSE(10)
107CLOSE(11)
108
109
110END PROGRAM shp_fxtr
111
Fortran 2003 interface to the shapelib http://shapelib.maptools.org/ library.
Definition shapelib.F90:53