convert coordinates to WGS84 using PROJ-library

src/Makefile

@@ -9,6 +9,7 @@ INCPATH = -I $(YAC_ROOT)/include \
 
 LIBS = \
         -lyac_utils \
+        -lproj \
         $(shell export PKG_CONFIG_PATH=$(YAC_ROOT)/lib/pkgconfig ; pkg-config yac --libs)
 
 CC = mpicc

src/elmer_grid.c

@@ -3,6 +3,8 @@
 #include <string.h>
 #include <mpi.h>
 
+#include <proj.h>
+
 #include "elmer_grid.h"
 
 struct glb2loc {
@@ -13,6 +15,31 @@ static int compare_glb2loc_glb (const void * a, const void * b) {
   return ((struct glb2loc*)a)->global_id - ((struct glb2loc*)b)->global_id;
 }
 
+static void convert2rad(
+  double * x_vertices, double * y_vertices, int nbr_vertices) {
+
+  // define transformation
+  PJ * P =
+    proj_create_crs_to_crs(
+      PJ_DEFAULT_CTX, "EPSG:3413", "+proj=longlat +datum=WGS84", NULL);
+
+  if (!P) {
+    fputs("failed to create transformation", stderr);
+    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
+  }
+
+  // transform all vertices
+  for (int i = 0; i < nbr_vertices; ++i) {
+    PJ_COORD src_coord = proj_coord(x_vertices[i], y_vertices[i], 0, 0);
+    PJ_COORD tgt_coord = proj_trans(P, PJ_FWD, src_coord);
+    x_vertices[i] = proj_torad(tgt_coord.lp.lam);
+    y_vertices[i] = proj_torad(tgt_coord.lp.phi);
+  }
+
+  // clean up
+  proj_destroy(P);
+}
+
 void read_grid(
   char const * grid_dir, int rank,
   int * nbr_vertices, int * nbr_cells, int ** num_vertices_per_cell,
@@ -21,7 +48,7 @@ void read_grid(
 
   if (!grid_dir) {
 
-    fputs("invalid grid_dir", stderr);
+    fputs("invalid grid_dir\n", stderr);
     MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
   }
 
@@ -67,6 +94,9 @@ void read_grid(
     glb2loc_vert[i].local_id = i;
   }
 
+  // convert coordiantes to radian
+  convert2rad(*x_vertices, *y_vertices, *nbr_vertices);
+
   enum {NUM_VERT_PER_CELL = 3};
   *num_vertices_per_cell =
     malloc((size_t)*nbr_cells * sizeof(**num_vertices_per_cell));
@@ -104,7 +134,7 @@ void read_grid(
     if ((j >= *nbr_vertices) ||
         (glb2loc_cell_vert[i].global_id != glb2loc_vert[j].global_id)) {
 
-      fputs("could not match cell vertices to list of vertices", stderr);
+      fputs("could not match cell vertices to list of vertices\n", stderr);
       MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
     }