1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 //
3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
4 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
5 // Use, modification and distribution is subject to the Boost Software License,
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
8 //
9 // SHAPELIB example
10
11 // Shapelib is a well-known and often used library to read (and write) shapefiles by Frank Warmerdam
12
13 // To build and run this example:
14 // 1) download shapelib from http://shapelib.maptools.org/
15 // 2) extract and put the source "shpopen.cpp" in project or makefile
16 // 3) download a shapefile, for example world countries from http://aprsworld.net/gisdata/world
17 // Alternativelly, install Shapelib using OSGeo4W installer from http://trac.osgeo.org/osgeo4w/
18 // that provides Windows binary packages
19
20 #include "shapefil.h"
21
22 #include <boost/geometry/geometry.hpp>
23 #include <boost/geometry/geometries/geometries.hpp>
24 #include <boost/geometry/geometries/point_xy.hpp>
25
26 using namespace boost::geometry;
27
28 template <typename T, typename F>
read_shapefile(const std::string & filename,std::vector<T> & polygons,F functor)29 void read_shapefile(const std::string& filename, std::vector<T>& polygons, F functor)
30 {
31 try
32 {
33 SHPHandle handle = SHPOpen(filename.c_str(), "rb");
34 if (handle <= 0)
35 {
36 throw std::string("File " + filename + " not found");
37 }
38
39 int nShapeType, nEntities;
40 double adfMinBound[4], adfMaxBound[4];
41 SHPGetInfo(handle, &nEntities, &nShapeType, adfMinBound, adfMaxBound );
42
43 for (int i = 0; i < nEntities; i++)
44 {
45 SHPObject* psShape = SHPReadObject(handle, i );
46
47 // Read only polygons, and only those without holes
48 if (psShape->nSHPType == SHPT_POLYGON && psShape->nParts == 1)
49 {
50 T polygon;
51 functor(psShape, polygon);
52 polygons.push_back(polygon);
53 }
54 SHPDestroyObject( psShape );
55 }
56 SHPClose(handle);
57 }
58 catch(const std::string& s)
59 {
60 throw s;
61 }
62 catch(...)
63 {
64 throw std::string("Other exception");
65 }
66 }
67
68
69 template <typename T>
convert(SHPObject * psShape,T & polygon)70 void convert(SHPObject* psShape, T& polygon)
71 {
72 double* x = psShape->padfX;
73 double* y = psShape->padfY;
74 for (int v = 0; v < psShape->nVertices; v++)
75 {
76 typename point_type<T>::type point;
77 assign_values(point, x[v], y[v]);
78 append(polygon, point);
79 }
80 }
81
82
main()83 int main()
84 {
85 std::string filename = "c:/data/spatial/shape/world_free/world.shp";
86
87 typedef model::polygon<model::d2::point_xy<double> > polygon_2d;
88 std::vector<polygon_2d> polygons;
89
90 try
91 {
92 read_shapefile(filename, polygons, convert<polygon_2d>);
93 }
94 catch(const std::string& s)
95 {
96 std::cout << s << std::endl;
97 return 1;
98 }
99
100 // Do something with the polygons, for example simplify them
101 for (std::vector<polygon_2d>::iterator it = polygons.begin(); it != polygons.end(); it++)
102 {
103 polygon_2d p;
104 simplify(*it, p, 0.01);
105 std::cout << it->outer().size() << "," << p.outer().size() << std::endl;
106 *it = p;
107 }
108 std::cout << "Simplified " << polygons.size() << std::endl;
109
110 double sum = 0;
111 for (std::vector<polygon_2d>::const_iterator it = polygons.begin(); it != polygons.end(); it++)
112 {
113 sum += area(*it);
114 }
115 std::cout << "Total area of " << polygons.size() << " polygons, total: " << sum << std::endl;
116
117 return 0;
118 }
119