I have a large number of objects defined by WGS84 coordinates which I would like to transform into Mercator space. There are other operations not shown in the MWE which prevent me from just doing everything in a single-threaded manner.
My current code is as follows:
//Compile with: g++ -O3 temp.cpp -lproj -fopenmp
#include <proj_api.h>
#include <cstdlib>
#include <vector>
#include <utility>
#include <cmath>
typedef std::vector<double> dvec;
typedef std::pair<dvec,dvec> dvpair;
//Define these here so we're not making them every time we run `pj_transform()`
const char* wgs84_str = "+init=epsg:4326";
const char* merc_str = "+init=epsg:3857";
const projPJ pj_wgs84 = pj_init_plus(wgs84_str);
const projPJ pj_merc = pj_init_plus(merc_str);
//Use the global projections to do projecty stuff
void ToMercator(std::vector<double> &x, std::vector<double> &y){
pj_transform(pj_wgs84, pj_merc, x.size(), 1, x.data(), y.data(), NULL);
}
int main(int argc, char **argv){
std::vector<dvpair> dvs(1000);
for(int i=0;i<dvs.size();i++){
dvs[i].first.push_back(2*M_PI*rand()-M_PI); //Longitude in [-180,180]
dvs[i].second.push_back(M_PI*rand()-M_PI/2); //Latitude in [-90,90]
}
#pragma omp parallel for
for(int i=0;i<dvs.size();i++)
ToMercator(dvs[i].first,dvs[i].second);
}
But it segfaults with a heisenbug. Getting a core dump with ulimit -c unlimited shows that pj_transform is raising the problem. This isn't too mysterious since some digging showed that proj4 is not thread safe.
The way to fix this is to create per-thread contexts for proj4 per this site. And this brings us to the meet of the problem.
Each thread needs to run pj_ctx_alloc() which will return a thread-specific variable of type projCtx, which then needs to be used with pj_init_plus_ctx() to make the projPJ objects responsible for the coordinate transforms.
How can I create thread-specific globals in OpenMP?
(You will want to answer that "globals are bad". I understand this, but note that removing them doesn't fundamentally change the question. I also have an aversion to "tramp data". In my actual application there are at least six layers of intermediate functions before I reach pj_transform(): I have a preference for not passing parameters that deep in this instance. You may disagree with the design choice, but this is not the right venue to argue about it.)
Edit This also did not work:
#include <proj_api.h>
#include <cstdlib>
#include <vector>
#include <utility>
#include <cmath>
typedef std::vector<double> dvec;
typedef std::pair<dvec,dvec> dvpair;
//Define these here so we're not making them every time we run `pj_transform()`
const char* wgs84_str = "+init=epsg:4326";
const char* merc_str = "+init=epsg:3857";
projPJ pj_wgs84;
projPJ pj_merc;
projCtx pj_ctx;
//Use the global projections to do projecty stuff
void ToMercator(std::vector<double> &x, std::vector<double> &y){
pj_transform(pj_wgs84, pj_merc, x.size(), 1, x.data(), y.data(), NULL);
}
int main(int argc, char **argv){
std::vector<dvpair> dvs(1000);
for(int i=0;i<dvs.size();i++){
dvs[i].first.push_back(2*M_PI*rand()-M_PI); //Longitude in [-180,180]
dvs[i].second.push_back(M_PI*rand()-M_PI/2); //Latitude in [-90,90]
}
#pragma omp parallel private(pj_ctx, pj_wgs84, pj_merc)
{
pj_ctx = pj_ctx_alloc();
pj_wgs84 = pj_init_plus_ctx(pj_ctx,wgs84_str);
pj_merc = pj_init_plus_ctx(pj_ctx,merc_str);
#pragma omp for
for(int i=0;i<dvs.size();i++)
ToMercator(dvs[i].first,dvs[i].second);
}
}
Edit
Using:
projPJ pj_wgs84;
projPJ pj_merc;
projCtx pj_ctx;
#pragma omp threadprivate(pj_wgs84)
#pragma omp threadprivate(pj_merc)
#pragma omp threadprivate(pj_ctx)
...
#pragma omp parallel private(pj_ctx, pj_wgs84, pj_merc)
{
pj_ctx = pj_ctx_alloc();
...
Gives a linker error:
/usr/bin/ld: pj_merc: TLS definition in /tmp/ccqZRl8N.o section .tbss mismatches non-TLS definition in /usr/lib/gcc/x86_64-linux-gnu/6/../../../x86_64-linux-gnu/libproj.so section .text
/usr/lib/gcc/x86_64-linux-gnu/6/../../../x86_64-linux-gnu/libproj.so: error adding symbols: Bad value
collect2: error: ld returned 1 exit status
