0

我有大量由 WGS84 坐标定义的对象，我想将它们转换为墨卡托空间。MWE 中未显示其他操作，这使我无法仅以单线程方式执行所有操作。

我目前的代码如下：

//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);
}

但它与heisenbug有段错误。获得带有ulimit -c unlimited显示pj_transform问题的核心转储。这并不太神秘，因为一些挖掘表明 proj4 不是线程安全的。

解决此问题的方法是为每个站点的 proj4 创建每个线程上下文。这使我们遇到了问题。

每个线程都需要运行pj_ctx_alloc()，这将返回一个特定于线程的变量 type projCtx，然后需要使用它pj_init_plus_ctx()来使projPJ对象负责坐标变换。

如何在 OpenMP 中创建特定于线程的全局变量？

（您会想回答“全局变量不好”。我理解这一点，但请注意，删除它们并不会从根本上改变问题。我也厌恶“流浪数据”。在我的实际应用中，至少有六个在我到达之前的中间函数层pj_transform()：在这种情况下，我倾向于不传递那么深的参数。您可能不同意设计选择，但这不是争论它的正确场所。）

编辑这也不起作用：

#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);
  }
}

编辑

使用：

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();

...

给出链接器错误：

/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
4

1 回答 1

0

以下作品：

#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";

//Using static means that these variables are file-scope. If we didn't use
//static then the variables might be accessed across many files and OpenMP would
//not like that
static projPJ  pj_wgs84;
static projPJ  pj_merc;
static projCtx pj_ctx;

//Set it so each thread has its own, private copy
#pragma omp threadprivate(pj_wgs84)
#pragma omp threadprivate(pj_merc)
#pragma omp threadprivate(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(10000);
  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]
  }

  //Do this somewhere convenient ... like here
  #pragma omp parallel
  {
    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 parallel for
  for(int i=0;i<dvs.size();i++)
    ToMercator(dvs[i].first,dvs[i].second);
}

悲伤的消息：在具有更多内核（24 对 4）的机器上测试它仍然会导致段错误。

于 2017-04-28T22:28:10.183 回答