/***************************************************************************** * * This file is part of Mapnik (c++ mapping toolkit) * * Copyright (C) 2014 Artem Pavlenko * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA * *****************************************************************************/ #include #include #include #include #include "quadtree.hpp" #include "shapefile.hpp" #include "shape_io.hpp" #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused-parameter" #pragma GCC diagnostic ignored "-Wunused-local-typedef" #include #include #include #pragma GCC diagnostic pop const int DEFAULT_DEPTH = 8; const double DEFAULT_RATIO=0.55; int main (int argc,char** argv) { using namespace mapnik; namespace po = boost::program_options; using std::string; using std::vector; using std::clog; using std::endl; bool verbose=false; unsigned int depth=DEFAULT_DEPTH; double ratio=DEFAULT_RATIO; vector shape_files; try { po::options_description desc("shapeindex utility"); desc.add_options() ("help,h", "produce usage message") ("version,V","print version string") ("verbose,v","verbose output") ("depth,d", po::value(), "max tree depth\n(default 8)") ("ratio,r",po::value(),"split ratio (default 0.55)") ("shape_files",po::value >(),"shape files to index: file1 file2 ...fileN") ; po::positional_options_description p; p.add("shape_files",-1); po::variables_map vm; po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm); po::notify(vm); if (vm.count("version")) { clog<<"version 0.3.0" <(); } if (vm.count("ratio")) { ratio = vm["ratio"].as(); } if (vm.count("shape_files")) { shape_files=vm["shape_files"].as< vector >(); } } catch (std::exception const& ex) { clog << "Error: " << ex.what() << endl; return -1; } clog << "max tree depth:" << depth << endl; clog << "split ratio:" << ratio << endl; //vector::const_iterator itr = shape_files.begin(); if (shape_files.size() == 0) { clog << "no shape files to index" << endl; return 0; } for (auto const& filename : shape_files) { clog << "processing " << filename << endl; std::string shapename (filename); boost::algorithm::ireplace_last(shapename,".shp",""); std::string shapename_full (shapename + ".shp"); if (! mapnik::util::exists (shapename_full)) { clog << "Error : file " << shapename_full << " does not exist" << endl; continue; } shape_file shp (shapename_full); if (! shp.is_open()) { clog << "Error : cannot open " << shapename_full << endl; continue; } int code = shp.read_xdr_integer(); //file_code == 9994 clog << code << endl; shp.skip(5*4); int file_length=shp.read_xdr_integer(); int version=shp.read_ndr_integer(); int shape_type=shp.read_ndr_integer(); box2d extent; shp.read_envelope(extent); clog << "length=" << file_length << endl; clog << "version=" << version << endl; clog << "type=" << shape_type << endl; clog << "extent:" << extent << endl; int pos=50; shp.seek(pos*2); quadtree tree(extent,depth,ratio); int count=0; while (true) { long offset=shp.pos(); int record_number=shp.read_xdr_integer(); int content_length=shp.read_xdr_integer(); shape_type = shp.read_ndr_integer(); box2d item_ext; if (shape_type==shape_io::shape_null) { if (pos >= file_length) { break; } else { // still need to increment pos, or the pos counter // won't indicate EOF until too late. pos+=4+content_length; continue; } } else if (shape_type==shape_io::shape_point) { double x=shp.read_double(); double y=shp.read_double(); item_ext=box2d(x,y,x,y); } else if (shape_type==shape_io::shape_pointm) { double x=shp.read_double(); double y=shp.read_double(); // skip m shp.read_double(); item_ext=box2d(x,y,x,y); } else if (shape_type==shape_io::shape_pointz) { double x=shp.read_double(); double y=shp.read_double(); // skip z shp.read_double(); // According to ESRI shapefile doc // A PointZ consists of a triplet of double-precision coordinates in the order X, Y, Z plus a // measure. // PointZ // { // Double X // X coordinate // Double Y // Y coordinate // Double Z // Z coordinate // Double M // Measure // } // But OGR creates shapefiles with M missing so we need to skip M only if present // NOTE: content_length is in 16-bit words if ( content_length == 18) { shp.read_double(); } item_ext=box2d(x,y,x,y); } else { shp.read_envelope(item_ext); shp.skip(2*content_length-4*8-4); } tree.insert(offset,item_ext); if (verbose) { clog << "record number " << record_number << " box=" << item_ext << endl; } pos+=4+content_length; ++count; if (pos >= file_length) break; } clog << " number shapes=" << count << endl; std::fstream file((shapename+".index").c_str(), std::ios::in | std::ios::out | std::ios::trunc | std::ios::binary); if (!file) { clog << "cannot open index file for writing file \"" << (shapename+".index") << "\"" << endl; } else { tree.trim(); std::clog<<" number nodes="<