rework projection code

- make use of osmium's proj wrapper class CRS
- use virtual functions instead of switching over projection
- remove duplicate code (mainly WSG84 to mercator transformation)
- use custom projections also with --proj 4326 and 900913 (fixes #514)
This commit is contained in:
Sarah Hoffmann
2016-01-12 00:01:07 +01:00
parent ccd4ebe9a6
commit 0a6238da51
18 changed files with 184 additions and 280 deletions

View File

@ -35,7 +35,7 @@ class GeometryFactory;
class CoordinateSequence;
}}
struct reprojection;
class reprojection;
class geometry_builder
{

View File

@ -14,7 +14,7 @@
std::shared_ptr<geometry_processor> geometry_processor::create(const std::string &type,
const options_t *options) {
std::shared_ptr<geometry_processor> ptr;
int srid = options->projection->project_getprojinfo()->srs;
int srid = options->projection->target_srs();
if (type == "point") {
ptr = std::make_shared<processor_point>(srid);

View File

@ -260,7 +260,7 @@ std::string database_options_t::conninfo() const
}
options_t::options_t():
prefix("planet_osm"), scale(DEFAULT_SCALE), projection(new reprojection(PROJ_SPHERE_MERC)), append(false), slim(false),
prefix("planet_osm"), scale(DEFAULT_SCALE), projection(reprojection::create_projection(PROJ_SPHERE_MERC)), append(false), slim(false),
cache(800), tblsmain_index(boost::none), tblsslim_index(boost::none), tblsmain_data(boost::none), tblsslim_data(boost::none), style(OSM2PGSQL_DATADIR "/default.style"),
expire_tiles_zoom(-1), expire_tiles_zoom_min(-1), expire_tiles_filename("dirty_tiles"), hstore_mode(HSTORE_NONE), enable_hstore_index(false),
enable_multi(false), hstore_columns(), keep_coastlines(false), parallel_indexing(true),
@ -319,13 +319,13 @@ options_t::options_t(int argc, char *argv[]): options_t()
keep_coastlines = true;
break;
case 'l':
projection.reset(new reprojection(PROJ_LATLONG));
projection.reset(reprojection::create_projection(PROJ_LATLONG));
break;
case 'm':
projection.reset(new reprojection(PROJ_SPHERE_MERC));
projection.reset(reprojection::create_projection(PROJ_SPHERE_MERC));
break;
case 'E':
projection.reset(new reprojection(-atoi(optarg)));
projection.reset(reprojection::create_projection(atoi(optarg)));
break;
case 'p':
prefix = optarg;
@ -492,7 +492,7 @@ options_t::options_t(int argc, char *argv[]): options_t()
//NOTE: this is hugely important if you set it inappropriately and are are caching nodes
//you could get overflow when working with larger coordinates (mercator) and larger scales
scale = (projection->get_proj_id() == PROJ_LATLONG) ? 10000000 : 100;
scale = (projection->target_latlon()) ? 10000000 : 100;
}
void options_t::check_options()

View File

@ -63,8 +63,8 @@ int main(int argc, char *argv[])
osmdata_t osmdata(middle, outputs);
fprintf(stderr, "Using projection SRS %d (%s)\n",
options.projection->project_getprojinfo()->srs,
options.projection->project_getprojinfo()->descr );
options.projection->target_srs(),
options.projection->target_desc());
//start it up
time_t overall_start = time(nullptr);

View File

@ -14,8 +14,6 @@
#include <iostream>
#include <memory>
#define SRID (reproj->project_getprojinfo()->srs)
#define CREATE_PLACE_TABLE \
"CREATE TABLE place (" \
" osm_type CHAR(1) NOT NULL," \
@ -608,10 +606,10 @@ int output_gazetteer_t::connect() {
int output_gazetteer_t::start()
{
reproj = m_options.projection;
int srid = m_options.projection->target_srs();
builder.set_exclude_broken_polygon(m_options.excludepoly);
places.srid_str = (boost::format("SRID=%1%;") % SRID).str();
places.srid_str = (boost::format("SRID=%1%;") % srid).str();
if(connect())
util::exit_nicely();
@ -638,7 +636,7 @@ int output_gazetteer_t::start()
pgsql_exec(Connection, PGRES_COMMAND_OK, CREATE_PLACE_ID_INDEX, "", "");
}
pgsql_exec(Connection, PGRES_TUPLES_OK, "SELECT AddGeometryColumn('place', 'geometry', %d, 'GEOMETRY', 2)", SRID);
pgsql_exec(Connection, PGRES_TUPLES_OK, "SELECT AddGeometryColumn('place', 'geometry', %d, 'GEOMETRY', 2)", srid);
pgsql_exec(Connection, PGRES_COMMAND_OK, "ALTER TABLE place ALTER COLUMN geometry SET NOT NULL");
}

View File

@ -47,8 +47,6 @@
#define BOOST_DIAGNOSTIC_INFO(e) boost::diagnostic_information((e))
#endif
#define SRID (reproj->project_getprojinfo()->srs)
/* example from: pg_dump -F p -t planet_osm gis
COPY planet_osm (osm_id, name, place, landuse, leisure, "natural", man_made, waterway, highway, railway, amenity, tourism, learning, building, bridge, layer, way) FROM stdin;
17959841 \N \N \N \N \N \N \N bus_stop \N \N \N \N \N \N -\N 0101000020E610000030CCA462B6C3D4BF92998C9B38E04940
@ -89,7 +87,7 @@ int output_pgsql_t::pgsql_out_way(osmid_t id, taglist_t &outtags,
{
/* Split long ways after around 1 degree or 100km */
double split_at;
if (m_options.projection->get_proj_id() == PROJ_LATLONG)
if (m_options.projection->target_latlon())
split_at = 1;
else
split_at = 100 * 1000;
@ -145,7 +143,7 @@ int output_pgsql_t::pgsql_out_relation(osmid_t id, const taglist_t &rel_tags,
}
/* Split long linear ways after around 1 degree or 100km (polygons not effected) */
if (m_options.projection->get_proj_id() == PROJ_LATLONG)
if (m_options.projection->target_latlon())
split_at = 1;
else
split_at = 100 * 1000;
@ -663,7 +661,8 @@ output_pgsql_t::output_pgsql_t(const middle_query_t* mid_, const options_t &opti
//have a different tablespace/hstores/etc per table
m_tables.push_back(std::shared_ptr<table_t>(
new table_t(
m_options.database_options.conninfo(), name, type, columns, m_options.hstore_columns, SRID,
m_options.database_options.conninfo(), name, type, columns, m_options.hstore_columns,
reproj->target_srs(),
m_options.append, m_options.slim, m_options.droptemp, m_options.hstore_mode,
m_options.enable_hstore_index, m_options.tblsmain_data, m_options.tblsmain_index
)

View File

@ -142,16 +142,13 @@ void parse_osmium_t::node(osmium::Node& node)
}
if (!m_bbox || m_bbox->contains(node.location())) {
double lat = node.location().lat_without_check();
double lon = node.location().lon_without_check();
m_proj->reproject(&lat, &lon);
auto c = m_proj->reproject(node.location());
convert_tags(node);
if (m_append) {
m_data->node_modify(node.id(), lat, lon, tags);
m_data->node_modify(node.id(), c.y, c.x, tags);
} else {
m_data->node_add(node.id(), lat, lon, tags);
m_data->node_add(node.id(), c.y, c.x, tags);
}
m_stats.add_node(node.id());
}

View File

@ -35,7 +35,7 @@
#include <osmium/handler.hpp>
struct reprojection;
class reprojection;
class osmdata_t;
class parse_stats_t

View File

@ -7,228 +7,130 @@
#include "config.h"
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <proj_api.h>
#include "reprojection.hpp"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/** must match expire.tiles.c */
#define EARTH_CIRCUMFERENCE 40075016.68
Projection_Info::Projection_Info(const char *descr_, const char *proj4text_, int srs_, const char *option_)
: descr(descr_), proj4text(proj4text_), srs(srs_), option(option_) {
}
namespace {
const struct Projection_Info Projection_Infos[] = {
/*PROJ_LATLONG*/ Projection_Info(
/*descr */ "Latlong",
/*proj4text*/ "+init=epsg:4326",
/*srs */ 4326,
/*option */ "-l" ),
/*PROJ_SPHERE_MERC*/ Projection_Info(
/*descr */ "Spherical Mercator",
/*proj4text*/ "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs",
/*srs */ 900913,
/*option */ "-m" )
void latlon2merc(double *lat, double *lon)
{
if (*lat > 85.07)
*lat = 85.07;
else if (*lat < -85.07)
*lat = -85.07;
using namespace osmium::geom;
*lon = (*lon) * EARTH_CIRCUMFERENCE / 360.0;
*lat = log(tan(PI/4.0 + deg_to_rad(*lat) / 2.0)) * EARTH_CIRCUMFERENCE/(PI*2);
}
class latlon_reprojection_t : public reprojection
{
public:
osmium::geom::Coordinates reproject(osmium::Location loc) const override
{
return osmium::geom::Coordinates(loc.lon_without_check(),
loc.lat_without_check());
}
void target_to_tile(double *lat, double *lon) const override
{
latlon2merc(lat, lon);
}
int target_srs() const override { return PROJ_LATLONG; }
const char *target_desc() const override { return "Latlong"; }
};
} // anonymous namespace
/* Positive numbers refer the to the table above, negative numbers are
assumed to refer to EPSG codes and it uses the proj4 to find those. */
reprojection::reprojection(int proj)
: Proj(proj), pj_source(nullptr), pj_target(nullptr), pj_tile(nullptr),
custom_projection(nullptr)
class merc_reprojection_t : public reprojection
{
char buffer[32];
/* hard-code the source projection to be lat/lon, since OSM XML always
* has coordinates in degrees. */
pj_source = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs");
/* hard-code the tile projection to be spherical mercator always.
* theoretically this could be made selectable but not all projections
* lend themselves well to making tiles; non-spherical mercator tiles
* are uncharted waters in OSM. */
pj_tile = pj_init_plus(Projection_Infos[PROJ_SPHERE_MERC].proj4text);
/* now set the target projection - the only one which is really variable */
if (proj >= 0 && proj < PROJ_COUNT)
{
pj_target = pj_init_plus(Projection_Infos[proj].proj4text);
}
else if (proj < 0)
{
if (snprintf(buffer, sizeof(buffer), "+init=epsg:%d", -proj ) >= (int)sizeof(buffer))
{
fprintf(stderr, "Buffer overflow computing proj4 initialisation string\n");
exit(1);
}
pj_target = pj_init_plus(buffer);
if (!pj_target)
{
fprintf (stderr, "Couldn't read EPSG definition (do you have /usr/share/proj/epsg?)\n");
exit(1);
}
}
if (!pj_source || !pj_target || !pj_tile)
{
fprintf(stderr, "Projection code failed to initialise\n");
exit(1);
}
if (proj >= 0)
return;
if (snprintf(buffer, sizeof(buffer), "EPSG:%d", -proj) >= (int)sizeof(buffer))
{
fprintf(stderr, "Buffer overflow computing projection description\n");
exit(1);
}
custom_projection = new Projection_Info(
strdup(buffer),
pj_get_def(pj_target, 0),
-proj, "-E");
}
reprojection::~reprojection()
{
pj_free(pj_source);
pj_source = nullptr;
pj_free(pj_target);
pj_target = nullptr;
pj_free(pj_tile);
pj_tile = nullptr;
if (custom_projection != nullptr) {
delete custom_projection;
}
}
struct Projection_Info const *reprojection::project_getprojinfo(void)
{
if( Proj >= 0 )
return &Projection_Infos[Proj];
else
return custom_projection;
}
void reprojection::reproject(double *lat, double *lon) const
{
double x[1], y[1], z[1];
/** Caution: This section is only correct if the source projection is lat/lon;
* so even if it looks like pj_source was just a variable, things break if
* pj_source is something else than lat/lon. */
if (Proj == PROJ_LATLONG)
return;
if (Proj == PROJ_SPHERE_MERC)
public:
osmium::geom::Coordinates reproject(osmium::Location loc) const override
{
/* The latitude co-ordinate is clipped at slightly larger than the 900913 'world'
* extent of +-85.0511 degrees to ensure that the points appear just outside the
* edge of the map. */
double lat = loc.lat_without_check();
double lon = loc.lon_without_check();
if (*lat > 85.07)
*lat = 85.07;
if (*lat < -85.07)
*lat = -85.07;
latlon2merc(&lat, &lon);
*lat = log(tan(M_PI/4.0 + (*lat) * DEG_TO_RAD / 2.0)) * EARTH_CIRCUMFERENCE/(M_PI*2);
*lon = (*lon) * EARTH_CIRCUMFERENCE / 360.0;
return;
return osmium::geom::Coordinates(lon, lat);
}
x[0] = *lon * DEG_TO_RAD;
y[0] = *lat * DEG_TO_RAD;
z[0] = 0;
void target_to_tile(double *, double *) const override
{ /* nothing */ }
/** end of "caution" section. */
int target_srs() const override { return PROJ_SPHERE_MERC; }
const char *target_desc() const override { return "Spherical Mercator"; }
};
pj_transform(pj_source, pj_target, 1, 1, x, y, z);
class generic_reprojection_t : public reprojection
{
public:
generic_reprojection_t(int srs)
: pj_target(srs), pj_source(PROJ_LATLONG), pj_tile("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")
{}
*lat = y[0];
*lon = x[0];
}
osmium::geom::Coordinates reproject(osmium::Location loc) const override
{
using namespace osmium::geom;
return transform(pj_source, pj_target,
Coordinates(deg_to_rad(loc.lon_without_check()),
deg_to_rad(loc.lat_without_check())));
}
/**
* Converts from (target) coordinates to tile coordinates.
*
* The zoom level for the coordinates is explicitly given in the
* variable map_width.
void target_to_tile(double *lat, double *lon) const override
{
auto c = transform(pj_target, pj_tile, osmium::geom::Coordinates(*lon, *lat));
*lon = c.x;
*lat = c.y;
}
int target_srs() const override { return PROJ_SPHERE_MERC; }
const char *target_desc() const override { return "Spherical Mercator"; }
private:
osmium::geom::CRS pj_target;
/** The projection of the source data. Always lat/lon (EPSG:4326). */
osmium::geom::CRS pj_source;
/** The projection used for tiles. Currently this is fixed to be Spherical
* Mercator. You will usually have tiles in the same projection as used
* for PostGIS, but it is theoretically possible to have your PostGIS data
* in, say, lat/lon but still create tiles in Spherical Mercator.
*/
void reprojection::coords_to_tile(double *tilex, double *tiley, double lon, double lat,
int map_width)
osmium::geom::CRS pj_tile;
};
} // anonymous namespace
reprojection *reprojection::create_projection(int srs)
{
double x[1], y[1], z[1];
x[0] = lon;
y[0] = lat;
z[0] = 0;
if (Proj == PROJ_LATLONG)
switch (srs)
{
x[0] *= DEG_TO_RAD;
y[0] *= DEG_TO_RAD;
case PROJ_LATLONG: return new latlon_reprojection_t();
case PROJ_SPHERE_MERC: return new merc_reprojection_t();
}
/* since pj_tile is always spherical merc, don't bother doing anything if
* destination proj is the same. */
if (Proj != PROJ_SPHERE_MERC)
{
pj_transform(pj_target, pj_tile, 1, 1, x, y, z);
/** FIXME: pj_transform could fail if coordinates are outside +/- 85 degrees latitude */
}
/* if ever pj_tile were allowed to be PROJ_LATLONG then results would have to
* be divided by DEG_TO_RAD here. */
*tilex = map_width * (0.5 + x[0] / EARTH_CIRCUMFERENCE);
*tiley = map_width * (0.5 - y[0] / EARTH_CIRCUMFERENCE);
return new generic_reprojection_t(srs);
}
/**
* Converts from (target) coordinates to coordinates in the tile projection (EPSG:3857)
*
* Do not confuse with coords_to_tile which does *not* calculate coordinates in the
* tile projection, but tile coordinates.
*/
void reprojection::target_to_tile(double *lat, double *lon) const
void reprojection::coords_to_tile(double *tilex, double *tiley,
double lon, double lat, int map_width)
{
if (Proj == PROJ_SPHERE_MERC)
return;
target_to_tile(&lat, &lon);
if (Proj == PROJ_LATLONG)
{
if (*lat > 85.07)
*lat = 85.07;
if (*lat < -85.07)
*lat = -85.07;
*lon = (*lon) * EARTH_CIRCUMFERENCE / 360.0;
*lat = log(tan(M_PI/4.0 + (*lat) * DEG_TO_RAD / 2.0)) * EARTH_CIRCUMFERENCE/(M_PI*2);
return;
}
double z = 0;
pj_transform(pj_target, pj_tile, 1, 1, lon, lat, &z);
}
int reprojection::get_proj_id() const
{
return Proj;
*tilex = map_width * (0.5 + lon / EARTH_CIRCUMFERENCE);
*tiley = map_width * (0.5 - lat / EARTH_CIRCUMFERENCE);
}

View File

@ -10,45 +10,53 @@
#include <boost/noncopyable.hpp>
struct Projection_Info {
Projection_Info(const char *descr_, const char *proj4text_, int srs_, const char *option_);
#include <osmium/geom/projection.hpp>
#include <osmium/osm/location.hpp>
const char *descr;
const char *proj4text;
const int srs;
const char *option;
};
enum Projection { PROJ_LATLONG = 4326, PROJ_SPHERE_MERC = 900913 };
enum Projection { PROJ_LATLONG = 0, PROJ_SPHERE_MERC, PROJ_COUNT };
struct reprojection : public boost::noncopyable
class reprojection : public boost::noncopyable
{
explicit reprojection(int proj);
~reprojection();
public:
virtual ~reprojection() = default;
struct Projection_Info const* project_getprojinfo(void);
void reproject(double *lat, double *lon) const;
void target_to_tile(double *lat, double *lon) const;
void coords_to_tile(double *tilex, double *tiley, double lon, double lat, int map_width);
int get_proj_id() const;
private:
int Proj;
/** The projection of the source data. Always lat/lon (EPSG:4326). */
void *pj_source;
/** The target projection (used in the PostGIS tables). Controlled by the -l/-M/-m/-E options. */
void *pj_target;
/** The projection used for tiles. Currently this is fixed to be Spherical
* Mercator. You will usually have tiles in the same projection as used
* for PostGIS, but it is theoretically possible to have your PostGIS data
* in, say, lat/lon but still create tiles in Spherical Mercator.
/**
* Reproject from the source projection lat/lon (EPSG:4326)
* to target projection.
*/
void *pj_tile;
virtual osmium::geom::Coordinates reproject(osmium::Location loc) const = 0;
struct Projection_Info *custom_projection;
/**
* Converts coordinates from target projection to tile projection (EPSG:3857)
*
* Do not confuse with coords_to_tile which does *not* calculate coordinates in the
* tile projection, but tile coordinates.
*/
virtual void target_to_tile(double *lat, double *lon) const = 0;
/**
* Converts from target coordinates to tile coordinates.
*
* The zoom level for the coordinates is explicitly given in the
* variable map_width.
*/
void coords_to_tile(double *tilex, double *tiley,
double lon, double lat, int map_width);
virtual int target_srs() const = 0;
virtual const char *target_desc() const = 0;
bool target_latlon() const
{
return target_srs() == PROJ_LATLONG;
}
/**
* Create a reprojection object with target srs `srs`.
*
* The target projection (used in the PostGIS tables).
* Controlled by the -l/-m/-E options.
*/
static reprojection *create_projection(int srs);
};
#endif

View File

@ -134,21 +134,21 @@ void test_outputs()
int get_random_proj(std::vector<std::string>& args)
{
int proj = rand() % (PROJ_COUNT + 1);
int proj = rand() % 3;
switch(proj)
{
case PROJ_LATLONG:
case PROJ_SPHERE_MERC:
args.push_back(reprojection(proj).project_getprojinfo()->option);
break;
default:
case 1:
args.push_back("-l");
return PROJ_LATLONG;
case 2:
args.push_back("-m");
return PROJ_SPHERE_MERC;
}
args.push_back("--proj");
//nice contiguous block of valid epsgs here randomly use one of those..
proj = (rand() % (2962 - 2308)) + 2308;
args.push_back((boost::format("%1%") % proj).str());
proj = -proj;
break;
}
return proj;
}
@ -206,7 +206,7 @@ void test_random_perms()
args.push_back("osm2pgsql");
//pick a projection
options.projection.reset(new reprojection(get_random_proj(args)));
options.projection.reset(reprojection::create_projection(get_random_proj(args)));
//pick a style file
std::string style = get_random_string(15);

View File

@ -39,7 +39,7 @@ int main(int argc, char *argv[]) {
options.num_procs = 1;
options.slim = true;
options.projection.reset(new reprojection(PROJ_LATLONG));
options.projection.reset(reprojection::create_projection(PROJ_LATLONG));
options.output_backend = "multi";
options.style = "tests/test_output_multi_line_trivial.style.json";

View File

@ -44,7 +44,7 @@ void check_output_poly_trivial(bool enable_multi, std::shared_ptr<pg::tempdb> db
options.slim = 1;
options.enable_multi = enable_multi;
options.projection.reset(new reprojection(PROJ_LATLONG));
options.projection.reset(reprojection::create_projection(PROJ_LATLONG));
options.output_backend = "multi";
options.style = "tests/test_output_multi_poly_trivial.style.json";

View File

@ -39,7 +39,7 @@ int main(int argc, char *argv[]) {
options.num_procs = 1;
options.slim = true;
options.projection.reset(new reprojection(PROJ_LATLONG));
options.projection.reset(reprojection::create_projection(PROJ_LATLONG));
options.output_backend = "multi";
options.style = "tests/test_output_multi_tags.json";

View File

@ -67,7 +67,7 @@ void test_area_base(bool latlon, bool reproj, double expect_area_poly, double ex
options.style = "default.style";
options.prefix = "osm2pgsql_test";
if (latlon) {
options.projection.reset(new reprojection(PROJ_LATLONG));
options.projection.reset(reprojection::create_projection(PROJ_LATLONG));
}
if (reproj) {
options.reproject_area = true;

View File

@ -125,8 +125,8 @@ void test_latlong() {
options.slim = true;
options.style = "default.style";
options.projection.reset(new reprojection(PROJ_LATLONG));
options.scale = (options.projection->get_proj_id() == PROJ_LATLONG) ? 10000000 : 100;
options.projection.reset(reprojection::create_projection(PROJ_LATLONG));
options.scale = (options.projection->target_latlon()) ? 10000000 : 100;
auto out_test = std::make_shared<output_pgsql_t>(mid_pgsql.get(), options);

View File

@ -101,7 +101,7 @@ int main() {
options_t options;
auto projection = std::make_shared<reprojection>(PROJ_SPHERE_MERC);
std::shared_ptr<reprojection> projection(reprojection::create_projection(PROJ_SPHERE_MERC));
options.projection = projection;
auto out_test = std::make_shared<test_output_t>(options);

View File

@ -103,7 +103,7 @@ int main(int argc, char *argv[]) {
std::string inputfile = "tests/test_multipolygon.osm";
options_t options;
auto projection = std::make_shared<reprojection>(PROJ_SPHERE_MERC);
std::shared_ptr<reprojection> projection(reprojection::create_projection(PROJ_SPHERE_MERC));
options.projection = projection;
auto out_test = std::make_shared<test_output_t>(options);