// Boost.Geometry (aka GGL, Generic Geometry Library)
//
// Copyright (c) 2007-2011 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2011 Bruno Lalande, Paris, France.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
// SHAPELIB example

// Shapelib is a well-known and often used library to read (and write) shapefiles by Frank Warmerdam

// To build and run this example:
// 1) download shapelib from http://shapelib.maptools.org/
// 2) extract and put the source "shpopen.cpp" in project or makefile
// 3) download a shapefile, for example world countries from http://aprsworld.net/gisdata/world
// Alternativelly, install Shapelib using OSGeo4W installer from http://trac.osgeo.org/osgeo4w/
// that provides Windows binary packages

#include "shapefil.h"

#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/domains/gis/io/wkt/stream_wkt.hpp>

using namespace boost::geometry;

template <typename T, typename F>
void read_shapefile(const std::string& filename, std::vector<T>& polygons, F functor)
{
    try
    {
        SHPHandle handle = SHPOpen(filename.c_str(), "rb");
        if (handle <= 0)
        {
            throw std::string("File " + filename + " not found");
        }

        int nShapeType, nEntities;
        double adfMinBound[4], adfMaxBound[4];
        SHPGetInfo(handle, &nEntities, &nShapeType, adfMinBound, adfMaxBound );

        for (int i = 0; i < nEntities; i++)
        {
            SHPObject* psShape = SHPReadObject(handle, i );

            // Read only polygons, and only those without holes
            if (psShape->nSHPType == SHPT_POLYGON && psShape->nParts == 1)
            {
                T polygon;
                functor(psShape, polygon);
                polygons.push_back(polygon);
            }
            SHPDestroyObject( psShape );
        }
        SHPClose(handle);
    }
    catch(const std::string& s)
    {
        throw s;
    }
    catch(...)
    {
        throw std::string("Other exception");
    }
}


template <typename T>
void convert(SHPObject* psShape, T& polygon)
{
    double* x = psShape->padfX;
    double* y = psShape->padfY;
    for (int v = 0; v < psShape->nVertices; v++)
    {
        typename point_type<T>::type point;
        assign_values(point, x[v], y[v]);
        append(polygon, point);
    }
}


int main()
{
    std::string filename = "c:/data/spatial/shape/world_free/world.shp";

    typedef model::polygon<model::d2::point_xy<double> > polygon_2d;
    std::vector<polygon_2d> polygons;

    try
    {
        read_shapefile(filename, polygons, convert<polygon_2d>);
    }
    catch(const std::string& s)
    {
        std::cout << s << std::endl;
        return 1;
    }

    // Do something with the polygons, for example simplify them
    for (std::vector<polygon_2d>::iterator it = polygons.begin(); it != polygons.end(); it++)
    {
        polygon_2d p;
        simplify(*it, p, 0.01);
        std::cout << it->outer().size() << "," << p.outer().size() << std::endl;
        *it = p;
    }
    std::cout << "Simplified " << polygons.size() << std::endl;

    double sum = 0;
    for (std::vector<polygon_2d>::const_iterator it = polygons.begin(); it != polygons.end(); it++)
    {
        sum += area(*it);
    }
    std::cout << "Total area of " << polygons.size() << " polygons, total: " << sum << std::endl;

    return 0;
}