User:Dbenbenn/procshape.cc

Here is a simple C++ program I use for converting Shapefile output to SVG format. Specifically, after using shpdump to decode the Shapefile, and proj to project the data, I use this program to select the shapes I want and make the output SVG-ready.

Suppose you have projected data in data, a list of names for the shapes in names, and a list of shape numbers you want (1-based) in numbers. Use

procshape data numbers names > foo.svg

Here's the code:

#include <fstream>
#include <limits>
#include <list>
#include <map>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cmath>

// The input numbers are xxx.yy
#define PREC 100
#define PRECDIGS 2

/*
Usage: procshape shapes numbers names
shapes is the file that shpdump outputs.
numbers is the shape numbers to select.
names is the name of each shape.
*/

// A county can have multiple shapes.  E.g. individual islands.
class shape {
  public:
    shape(long _shapenum, long _parts, long _points);
    void add_point(long x, long y);
    void new_part(void);
    void optimize(void);
	// drop redundant points
    void output_shape(std::ostream &output) const;
    long get_shapenum(void) const { return shapenum; };

  private:
    std::list<std::list<std::pair<long, long> > > parts;
    std::list<std::list<std::pair<long, long> > >::iterator this_part;
    long shapenum;
    long expected_parts;
    long expected_points;
    long num_parts;
    long num_points;
};

shape::shape(long _shapenum, long _parts, long _points) :
	shapenum(_shapenum),
	expected_parts(_parts),
	expected_points(_points),
	num_parts(0),
	num_points(0)
{
    this_part = parts.end();
    new_part();
}

void shape::add_point(long x, long y)
{
    this_part->push_back(std::pair<long, long>(x, y));
    ++num_points;
}

void shape::new_part(void)
{
    parts.push_back(std::list<std::pair<long, long> >());
    ++this_part;
    ++num_parts;
}

void shape::optimize(void)
{
    for (std::list<std::list<std::pair<long, long> > >::iterator i = parts.begin(); i != parts.end(); ++i) {
	if (i->empty() || ++(i->begin()) == i->end())
	    continue;

	std::list<std::pair<long, long> >::const_iterator j1 = i->begin();
	std::list<std::pair<long, long> >::iterator j2 = ++(i->begin());
	std::list<std::pair<long, long> >::iterator j3 = j2;
	++j3;

	while (j3 != i->end()) {
	    if ((j2->second - j1->second) * (j3->first - j2->first) == (j3->second - j2->second) * (j2->first - j1->first))
		i->erase(j2++);
	    else {
		++j1;
		++j2;
	    }
	    ++j3;
	}
	j3 = i->begin();
	if ((j2->second - j1->second) * (j3->first - j2->first) == (j3->second - j2->second) * (j2->first - j1->first))
	    i->erase(j2);
	else
	    j1 = j2;
	j2 = j3;
	++j3;
	if ((j2->second - j1->second) * (j3->first - j2->first) == (j3->second - j2->second) * (j2->first - j1->first))
	    i->erase(j2);
    }
}

inline std::ostream& operator<<(std::ostream &out, const shape& shp)
{
    shp.output_shape(out);
    return out;
}

void shape::output_shape(std::ostream &out) const
{
    if (num_points != expected_points || num_parts != expected_parts)
	throw(42);

    std::list<std::list<std::pair<long, long> > >::const_iterator part;

    for (part = parts.begin(); part != parts.end(); ++part) {
	if (part->empty() || ++(part->begin()) == part->end())
	    throw(2);
	std::list<std::pair<long, long> >::const_iterator point = part->begin();

	out << "M " << point->first << ',' << point->second << " l";

	long lastx = point->first;
	long lasty = point->second;

	++point;
	for (; point != part->end(); ++point) {
	    out << ' ' << point->first - lastx << ',' << point->second - lasty;
	    lastx = point->first;
	    lasty = point->second;
	}
	out << " z";
    }
}

class County {
  public:
	// Skip ahead to shape number shapenum in input, and read it.
    void read_shape(std::istream &input, long shapenum);
    void write_county(std::ostream &output, std::string name) const;
    County(void);

	// Bounding box of all counties
    static long get_x(void) { return minx; };
    static long get_y(void) { return miny; };
    static long get_width(void) { return maxx - minx; };
    static long get_height(void) { return maxy - miny; };

  private:
    std::list<shape> shapes;
    long numshapes;

	// Bounding box for all counties that were stored (not skipped)
    static long minx;
    static long miny;
    static long maxx;
    static long maxy;

    long read_header(std::istream &input, long &numparts, long &numpoints);
};

long County::minx = std::numeric_limits<long>::max();
long County::miny = std::numeric_limits<long>::max();
long County::maxx = std::numeric_limits<long>::min();
long County::maxy = std::numeric_limits<long>::min();

County::County(void) :
	numshapes(0)
{
}

long County::read_header(std::istream &input, long &numparts, long &numpoints)
{
    std::string head;
    std::getline(input, head);
    std::istringstream format(head);
    std::string polygon, parts, points;
    long shapenum;

    format >> polygon >> shapenum >> parts >> numparts >> points >> numpoints;
    if (polygon.compare("polygon") != 0 || parts.compare("parts") != 0 ||
	    points.compare("points") != 0)
	throw(1);
    return shapenum;
}

void paranoid_check(const std::string &line, long x, long y)
{
    long absx = std::abs(x), absy = std::abs(y);
    std::ostringstream foo;
    if (x < 0 || x == 0 && line[0] == '-')
	foo << '-';
    foo << absx / PREC << '.' << std::setfill('0') << std::setw(PRECDIGS) << absx % PREC << std::setw(0)
	<< '\t';
    if (y < 0 || y == 0 && line[line.find('\t') + 1] == '-')
	foo << '-';
    foo << absy / PREC << '.' << std::setfill('0') << std::setw(PRECDIGS) << std::abs(absy) % PREC << std::setw(0);
    if (line.compare(foo.str()) != 0) {
	std::cerr << "Paranoid, line is " << line << ", foo is " << foo.str() << ", x is " << x << ", y is " << y << '\n';
//	exit(2);
    }
}

void County::read_shape(std::istream &input, long shapenum)
{
    long numparts, numpoints;

    while (shapenum != read_header(input, numparts, numpoints)) {
	for (long i = 0; i < numparts - 1 + numpoints; ++i) {
	    std::string line;
	    std::getline(input, line);
	}
    }
    numshapes++;
    shape newshape(shapenum, numparts, numpoints);

    for (long i = 0; i < numparts + numpoints - 1; ++i) {
	std::string line;
	std::getline(input, line);
	if (line.compare("part") == 0)
	    newshape.new_part();
	else {
	    double xdbl, ydbl;
	    long x, y;
	    std::istringstream foo(line);

	    foo >> xdbl >> ydbl;
	    x = (long) round(xdbl * PREC);
	    y = (long) round(-ydbl * PREC);

	    // These two lines are not equivalent to the above.  Consider 232.73
	    // x = (long) (PREC * xdbl);
	    // y = (long) (PREC * ydbl);
	    paranoid_check(line, x, -y);

	    minx = std::min(x, minx);
	    miny = std::min(y, miny);
	    maxx = std::max(x, maxx);
	    maxy = std::max(y, maxy);

	    newshape.add_point(x, y);
	}
    }
    newshape.optimize();
    shapes.push_back(newshape);
}

void County::write_county(std::ostream &output, std::string name) const
{
    if (numshapes == 0)
	throw(4);
    if (numshapes == 1) {
	output << "<path id=\"" << name << "\" d=\"" << *shapes.begin() << "\" />\n";
    } else {
#ifndef STATE
	output << "<g id=\"" << name << "\">\n";
#endif
	for (std::list<shape>::const_iterator i = shapes.begin(); i != shapes.end(); ++i)
	    output << "<path id=\"" << i->get_shapenum() << "\" d=\"" << *i << "\" />\n";
#ifndef STATE
	output << "</g>\n";
#endif
    }
}

void write_output_file(const std::map<std::string, County>& counties, std::ostream &output)
{
    long width = County::get_width();
    long height = County::get_height();
    long scalewidth = width;
    long widthremainder = 0;
    long scaleheight = height;
    long heightremainder = 0;
    long decpos = 1;
    int digs = 0;

    while (scalewidth > 12800 && scaleheight > 10240) {
	widthremainder += decpos * (scalewidth % 10);
	scalewidth /= 10;
	heightremainder += decpos * (scaleheight % 10);
	scaleheight /= 10;
	decpos *= 10;
	++digs;
    }

    long thick_line = (long) (0.007 * width);
    long thin_line = (long) (0.0016 * width);

#ifndef STATE
    output << "<?xml version=\"1.0\" standalone=\"yes\"?>\n"
	"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n"
	"<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" version=\"1.1\" width=\""
	<< scalewidth << "." << std::setfill('0')
	<< std::setw(digs) << widthremainder << std::setw(0)
	<< "\" height=\"" << scaleheight << "."
	<< std::setw(digs) << heightremainder << std::setw(0)
	<< "\" viewBox=\"" << County::get_x() << ',' << County::get_y() << ' '
	<< width << ',' << height << "\">\n"
	<< "<clipPath id=\"state_clip_path\">\n"
	"<g id=\"state_outline\">\n"
	"STATEHERE\n"
	"</g>\n"
	"</clipPath>\n"
	"<g stroke=\"black\" fill=\"none\" stroke-linejoin=\"round\" stroke-width=\""
	<< thin_line
	<< "\" clip-path=\"url(#state_clip_path)\">\n"
	"<use xlink:href=\"#state_outline\" fill=\"white\" stroke-width=\""
	<< thick_line
	<< "\" />\n";
#endif
    for (std::map<std::string, County>::const_iterator i = counties.begin(); i != counties.end(); ++i)
	i->second.write_county(output, i->first);
#ifndef STATE
    output << "</g>\n"
	"<use xlink:href=\"#FOO\" stroke=\"none\" fill=\"red\" />\n"
	"</svg>\n";
#endif
}

int main(int argc, char **argv)
{
    std::ifstream shapes_f(argv[1], std::ios::in);
    std::ifstream nums_f(argv[2], std::ios::in);
    std::ifstream names_f(argv[3], std::ios::in);
    long names_line = 0;
    std::map<std::string, County> counties;
    std::string line;

    std::getline(shapes_f, line);

    if (line.compare("type polygon") != 0) {
	std::cerr << "I can only deal with polygons right now\n";
	throw (1);
    }
    while (!nums_f.eof()) {
	long shapenum;
	std::string name;

	nums_f >> shapenum;
	if (nums_f.eof())
	    break;
	while (shapenum > names_line) {
	    std::getline(names_f, name);
	    ++names_line;
	}
	if (name.empty())
	    std::cerr << "Ignoring unnamed shape " << shapenum << '\n';
	else
	    counties[name].read_shape(shapes_f, shapenum);
    }
    write_output_file(counties, std::cout);

    return 0;
}