void coordconverter::convert(const std::vector<latlong>& latlongs, std::vector<coord>& coords)
{
  projPJ pjMerc = pj_init_plus("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs");
  if(pjMerc == NULL) throw hdsrvexception("pj_init_plus merc failed"); 

  projPJ pjLatLon = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs");
  if(pjLatLon == NULL) 
	{
		pj_free(pjMerc);
		throw hdsrvexception("pj_init_plus proj failed");
	}

	  
  for(std::vector<latlong>::const_iterator it = latlongs.begin(); it != latlongs.end(); it++)
  {
    double x = it->getLonDouble() * DEG_TO_RAD;
    double y = it->getLatDouble() * DEG_TO_RAD;

    if(pj_transform(pjLatLon, pjMerc, 1, 1, &x, &y, NULL) != 0)
    {
			pj_free(pjMerc);
			pj_free(pjLatLon);
			throw hdsrvexception("pj_transform failed");
    }
    coord c;
    c.setNDouble(y);
    c.setEDouble(x);
    coords.push_back(c);
  }
  
  pj_free(pjMerc);
  pj_free(pjLatLon);

    
}

int pj_init_and_transform( const char *from_projection, const char *to_projection, \
   long point_count, double *x, double *y) {

	 projPJ pj_from, pj_to;
	 int i;
	 int p;

    if (!(pj_from = pj_init_plus(from_projection)) ) {
	    printf("\nPROJ4 ERROR: There was a problem creating a PROJ4 object; "
	      "something is\nwrong with the PROJ4 string you supplied."
	      "\nOffending string: '%s'",from_projection);
        exit(1);
	}

    if (!(pj_to = pj_init_plus(to_projection)) ) {
	    printf("\nPROJ4 ERROR: There was a problem creating a PROJ4 object; "
	      "something is\nwrong with the PROJ4 string you supplied."
	      "\nOffending string: '%s'",to_projection);
        exit(1);
	}

    p = pj_transform(pj_from, pj_to, point_count, 1, x, y, NULL );

    pj_free(pj_to);
    pj_free(pj_from);

	return(p);
}

GeoConvHelper&
GeoConvHelper::operator=(const GeoConvHelper& orig) {
    myProjString = orig.myProjString;
    myOffset = orig.myOffset;
    myProjectionMethod = orig.myProjectionMethod;
    myOrigBoundary = orig.myOrigBoundary;
    myConvBoundary = orig.myConvBoundary;
    myGeoScale = orig.myGeoScale;
    myUseInverseProjection = orig.myUseInverseProjection;
#ifdef HAVE_PROJ
    if (myProjection != 0) {
        pj_free(myProjection);
        myProjection = 0;
    }
    if (myInverseProjection != 0) {
        pj_free(myInverseProjection);
        myInverseProjection = 0;
    }
    if (myGeoProjection != 0) {
        pj_free(myGeoProjection);
        myGeoProjection = 0;
    }
    if (orig.myProjection != 0) {
        myProjection = pj_init_plus(orig.myProjString.c_str());
    }
    if (orig.myInverseProjection != 0) {
        myInverseProjection = pj_init_plus(pj_get_def(orig.myInverseProjection, 0));
    }
    if (orig.myGeoProjection != 0) {
        myGeoProjection = pj_init_plus(pj_get_def(orig.myGeoProjection, 0));
    }
#endif
    return *this;
}

void GridGeometryTest::setUp()
{
	grid = new GridGeometry(
		"+proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +a=6367470.0 +no_defs",
		GridGeometry::LeftLowerHorizontal,
		248,
		400,
		0.1,// * DEG_TO_RAD,
		0.1,// * DEG_TO_RAD,
		5.75,// * DEG_TO_RAD,
		-13.25// * DEG_TO_RAD
	); // Hirlam 10 grid
	hirlam10Proj = pj_init_plus( "+proj=ob_tran +o_proj=longlat +lon_0=-40 +o_lat_p=22 +a=6367470.0 +no_defs" );
	if ( !hirlam10Proj )
		throw std::runtime_error( "Invalid PROJ definition for Hirlam10" );
	hirlam20Proj = pj_init_plus( "+proj=ob_tran +o_proj=longlat +lon_0=0 +o_lat_p=25 +a=6367470.0 +no_defs" );
	if ( !hirlam20Proj )
		throw std::runtime_error( "Invalid PROJ definition for Hirlam20" );
	targetProj = pj_init_plus( "+proj=longlat +ellps=WGS84 +no_defs" );
	if ( !targetProj )
		throw std::runtime_error( "Invalid PROJ definition for target" );
	utmProj = pj_init_plus( "+proj=utm +lon_0=15e +datum=WGS84 +units=m +no_defs" );
	if ( !utmProj )
		throw std::runtime_error( "Invalid PROJ definition for utm" );
}

static VALUE method_proj4_initialize_copy(VALUE self, VALUE orig)
{
    RGeo_Proj4Data* self_data;
    projPJ pj;
    RGeo_Proj4Data* orig_data;
    char* str;

    // Clear out any existing value
    self_data = RGEO_PROJ4_DATA_PTR(self);
    pj = self_data->pj;
    if (pj) {
        pj_free(pj);
        self_data->pj = NULL;
        self_data->original_str = Qnil;
    }

    // Copy value from orig
    orig_data = RGEO_PROJ4_DATA_PTR(orig);
    if (!NIL_P(orig_data->original_str)) {
        self_data->pj = pj_init_plus(RSTRING_PTR(orig_data->original_str));
    }
    else {
        str = pj_get_def(orig_data->pj, 0);
        self_data->pj = pj_init_plus(str);
        pj_dalloc(str);
    }
    self_data->original_str = orig_data->original_str;
    self_data->uses_radians = orig_data->uses_radians;

    return self;
}

/*!
 * \brief
 * executes reprojection 
 * 
 * JNI informations:
 * Class:     org_proj4_Projections
 * Method:    transform
 * Signature: ([D[D[DLjava/lang/String;Ljava/lang/String;JI)V
 * 
 *
 * \param env - parameter used by jni (see JNI specification)
 * \param parent - parameter used by jni (see JNI specification)
 * \param firstcoord - array of x coordinates
 * \param secondcoord - array of y coordinates
 * \param values - array of z coordinates
 * \param src - definition of the source projection
 * \param dest - definition of the destination projection
 * \param pcount
 * \param poffset
*/
JNIEXPORT void JNICALL Java_org_proj4_Projections_transform
  (JNIEnv * env, jobject parent, jdoubleArray firstcoord, jdoubleArray secondcoord, jdoubleArray values, jstring src, jstring dest, jlong pcount, jint poffset)
{
	int i;
	projPJ src_pj, dst_pj;
	char * srcproj_def = (char *) (*env)->GetStringUTFChars (env, src, 0); 
	char * destproj_def = (char *) (*env)->GetStringUTFChars (env, dest, 0);

	if (!(src_pj = pj_init_plus(srcproj_def)))
		exit(1);
	if (!(dst_pj = pj_init_plus(destproj_def)))
		exit(1);

	double *xcoord = (* env)-> GetDoubleArrayElements(env, firstcoord, NULL); 
	double *ycoord = (* env) -> GetDoubleArrayElements(env, secondcoord, NULL); 
	double *zcoord = (* env) -> GetDoubleArrayElements(env, values, NULL); 

        pj_transform( src_pj, dst_pj, pcount,poffset, xcoord, ycoord, zcoord);

	(* env)->ReleaseDoubleArrayElements(env,firstcoord,(jdouble *) xcoord,JNI_COMMIT);
	(* env)->ReleaseDoubleArrayElements(env,secondcoord,(jdouble *) ycoord,JNI_COMMIT);
	(* env)->ReleaseDoubleArrayElements(env,values,(jdouble *) zcoord,JNI_COMMIT);

	pj_free( src_pj );
	pj_free( dst_pj );
}

void Converter::setUTMZone(UtmRegion zone){
    const char* src_proj_str;
    const char* dst_proj_str;

    switch (zone) {
        case UtmRegion::Beijing:
            src_proj_str = BJ_SRC_PROJ;
            dst_proj_str = BJ_DST_PROJ;        
            break;
        case UtmRegion::California:
            src_proj_str = SF_SRC_PROJ;
            dst_proj_str = SF_DST_PROJ;
            break;
        default:
            break;
    }

    if (!(src_proj_ = pj_init_plus(SF_SRC_PROJ))) {
        fprintf(stderr, "Error! Cannot initialize latlon to XY projector!\n");
        exit(1);
    }
    if (!(dst_proj_ = pj_init_plus(SF_DST_PROJ))) {
        fprintf(stderr, "Error! Cannot initialize latlon to XY projector!\n");
        exit(1);
    }
}

// reproject shape
shape_t*
shape_proj(
    const shape_t* shape,
    const char*    from,
    const char*    to   ){

    shape_t* projected = shape_copy(shape);
    projPJ old_prj = pj_init_plus(from);
    projPJ new_prj = pj_init_plus(to);

    for(uint32_t i=0; i<kv_size(projected->points); i++) {
        point_t* p = &kv_A(projected->points, i);
        p->x *= DEG_TO_RAD;
        p->y *= DEG_TO_RAD;
        int32_t err = pj_transform(old_prj, new_prj, 1, 0, &p->x, &p->y, NULL);
        if (err)
            fprintf(stderr, "ERR%d %s\n", err, pj_strerrno(err));
        assert(err == 0);
    }

    for(uint32_t i=0; i<kv_size(projected->hull); i++) {
        point_t* p = &kv_A(projected->hull, i);
        p->x *= DEG_TO_RAD;
        p->y *= DEG_TO_RAD;
        int32_t err = pj_transform(old_prj, new_prj, 1, 0, &p->x, &p->y, NULL);
        if (err)
            fprintf(stderr, "ERR%d %s\n", err, pj_strerrno(err));
        assert(err == 0);
    }

    pj_free(old_prj);
    pj_free(new_prj);

    return projected;
}

char *msudf_transform(UDF_INIT *initid, UDF_ARGS *args, char *buf,
	unsigned long *length, char *is_null, char *error)
{
	unsigned char *wkb = NULL;
	size_t wkbsize;
	GEOSGeom g1,g2 =NULL;
	
	msudf_params *params;
	int srid_src,srid_dst;


	params = (msudf_params *) initid->ptr;
	

	DEBUG("msudf_transform");
	wkb = (unsigned char *) (args->args[0]);
	g1 = GEOSGeomFromWKB_buf(wkb + 4,args->lengths[0] - 4);
	srid_src = msudf_getInt(wkb);
	wkb = NULL;

	DEBUG("srid_src: %d",srid_src);
	if ((params->pj_src == NULL) || (params->srid_src != srid_src)) {
		params->pj_src = pj_init_plus(args->args[1]);
		params->srid_src = srid_src;
	}

	srid_dst = msudf_getInt((unsigned char *) args->args[2]);
	DEBUG("srid_dst: %d",srid_dst);
	if ((params->pj_dst == NULL) || (params->srid_dst != srid_dst)) {
		params->pj_dst = pj_init_plus(args->args[3]);
		params->srid_dst = srid_dst;
	}

	if (params->pj_src != NULL && params->pj_dst != NULL && g1 != NULL) {
		g2 = gu_transformGeom(g1,params->pj_src,params->pj_dst);
		wkb = GEOSGeomToWKB_buf(g2,&wkbsize);
	} else {
		// initid->ptr = NULL;
		strcpy(error,"Invalid geometry.");
	}

	if (g1 != NULL) GEOSGeom_destroy(g1);
	if (g2 != NULL) GEOSGeom_destroy(g2);


	if (wkb != NULL) {
		*length = (long)wkbsize + 4;
		if (params->buf != NULL) free(params->buf);
		params->buf = (char *) malloc(*length);
		memcpy(params->buf,args->args[2],4);
		memcpy((char *)params->buf + 4,wkb,wkbsize);
		GEOSFree((char *)wkb);
		return params->buf;
	} else {
		*is_null = 1;
		return NULL;
	}

}

ProjectionSettings::ProjectionSettings()
    : ui(new Ui::ProjectionSettings)
{
    inst = this;
    ui->setupUi(this);

    connect(this, SIGNAL(accepted()), this, SLOT(okPressed()));

    ui->XOffsetSpin->setDecimals(10);
    ui->XOffsetSpin->setMaximum(10000000);
    ui->XOffsetSpin->setMinimum(-10000000);
    //ui->XOffsetSpin->setValue(-5439122.807299255);
    ui->XOffsetSpin->setValue(-3506000);
    //ui->XOffsetSpin->setValue(-926151);
    ui->YOffsetSpin->setDecimals(10);
    ui->YOffsetSpin->setMaximum(10000000);
    ui->YOffsetSpin->setMinimum(-10000000);
    //ui->YOffsetSpin->setValue(-984970.1841083583);
    //ui->YOffsetSpin->setValue(-3463995);
    ui->YOffsetSpin->setValue(-5400147);
    ui->ZOffsetSpin->setDecimals(10);
    ui->ZOffsetSpin->setMaximum(10000000);
    ui->ZOffsetSpin->setMinimum(-10000000);
    //ui->ZOffsetSpin->setValue(-399.4944465);
    ui->ZOffsetSpin->setValue(0.0);
    ui->FromDatumEdit->setText(QString("WGS84"));
    ui->ToProjectionEdit->setText(QString("tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam"));
    //ui->ToProjectionEdit->setText(QString("tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +nadgrids=/data/porsche/BETA2007.gsb"));
    //ui->ToProjectionEdit->setText(QString("utm +zone=50 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"));
    QString projFromString = "+proj=latlong +datum=" + ui->FromDatumEdit->text();
    //std::string projToString = "+proj=merc +x_0=-1008832.89 +y_0=-6179385.47";.
    QString projToString = "+proj=" + ui->ToProjectionEdit->text();
    XOffset = ui->XOffsetSpin->value();
    YOffset = ui->YOffsetSpin->value();
    ZOffset = ui->ZOffsetSpin->value();
    
    projPJ new_pj_from, new_pj_to;
    if (!(new_pj_from = pj_init_plus(projFromString.toUtf8().constData())))
    {
        //std::cerr << "RoadSystem::parseIntermapRoad(): couldn't initalize projection source: " << projFromString << std::endl;
        //return false;
    }
    else
    {
        pj_from = new_pj_from;
    }
    if (!(new_pj_to = pj_init_plus(projToString.toUtf8().constData())))
    {
        //std::cerr << "RoadSystem::parseIntermapRoad(): couldn't initalize projection target: " << projToString << std::endl;
        //return false;
    }
    else
    {
        pj_to = new_pj_to;
    }
    inst = this;
}

// reproject shapes
shapes_v
shape_proj(
    const shapes_v* shapes, 
    const char* from, 
    const char* to){

    projPJ old_prj = pj_init_plus(from);
    projPJ new_prj = pj_init_plus(to);

    shapes_v shapes_prj;
    kv_init(shapes_prj);
    
    shapes_prj.min = (point_t){ INFINITY, INFINITY};
    shapes_prj.max = (point_t){-INFINITY,-INFINITY};
    double k = 0.0;
    for(uint32_t s=0; s<shapes->n; s++) {
        shape_v* shape = &shapes->a[s];
        shape_v shape_prj;
        kv_init(shape_prj);
        for(uint32_t p=0; p<shape->n; p++) {
            point_t pnt = shape->a[p];
            pnt.x *= DEG_TO_RAD;
            pnt.y *= DEG_TO_RAD;
            int32_t err = pj_transform(old_prj, new_prj, 1, 0, &pnt.x, &pnt.y, NULL);
            if (err) {
                fprintf(stderr, "ERR%d %s\n", err, pj_strerrno(err));
                continue;
            }

            // cumulitive average for center
            if(k>=1.0) {
                shapes_prj.center.x = (k-1.0)/k*shapes_prj.center.x + pnt.x/k;
                shapes_prj.center.y = (k-1.0)/k*shapes_prj.center.y + pnt.y/k;
            }else {
                shapes_prj.center.x = pnt.x;
                shapes_prj.center.y = pnt.y;
            }
            k+=1.0;

            // new bounds
            if(pnt.x>shapes_prj.max.x) shapes_prj.max.x = pnt.x;
            if(pnt.y>shapes_prj.max.y) shapes_prj.max.y = pnt.y;
            if(pnt.x<shapes_prj.min.x) shapes_prj.min.x = pnt.x;
            if(pnt.y<shapes_prj.min.y) shapes_prj.min.y = pnt.y;

            kv_push(point_t, shape_prj, pnt);
        }
        kv_push(shape_v, shapes_prj, shape_prj);
    }

    pj_free(old_prj);
    pj_free(new_prj);

    return shapes_prj;
}

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

static int test_nad27_ll(double lat, double lon,
                         double expected_lat, double expected_lon)
{
    double lats[1];
    double lons[1];
    double hts[1];

    char input_projection_info[255];
    char output_projection_info[255];

    projPJ input_projection, output_projection;

    sprintf(input_projection_info, "+proj=latlong +datum=NAD27");
    sprintf(output_projection_info, "+proj=latlong +datum=NAD83");

    input_projection = pj_init_plus(input_projection_info);
    output_projection = pj_init_plus(output_projection_info);
    
    lats[0] = lat * D2R;
    lons[0] = lon * D2R;
    hts[0] = 0;

    pj_transform (input_projection, output_projection, 1, 1, 
                  lats, lons, hts);

    if (pj_errno != 0) {
        asfPrintWarning("libproj error: %s\n", pj_strerrno(pj_errno));
    }

    pj_free(input_projection);
    pj_free(output_projection);

    lats[0] *= R2D;
    lons[0] *= R2D;

    CU_ASSERT(double_equals(lats[0], expected_lat, 6));
    CU_ASSERT(double_equals(lons[0], expected_lon, 6));

    if (double_equals(lats[0], expected_lat, 6) &&
        double_equals(lons[0], expected_lon, 6))
    {
        //asfPrintStatus("Proj (fwd): %f, %f ... OK\n", lat, lon);
        return TRUE;
    }
    else
    {
        asfPrintStatus("Proj (fwd): %f, %f ... ERROR\n", lat, lon);
        asfPrintStatus("Result:                 %.10f %.10f (%.10f)\n",
                       lats[0], lons[0], hts[0]);
        asfPrintStatus("Expected:               %.10f %.10f\n",
                       expected_lat, expected_lon);
        return FALSE;
    }
}

void process_projection(const projection_t *prj)
{

    projPJ pj_tgt, pj_latlong;
    const point_t *test_point;
    int i;
    double x,y;
    double x1,y1;

    pj_tgt=pj_init_plus(prj->init_string);
    pj_latlong=pj_init_plus("+proj=latlong +ellps=WGS84");

    if (!pj_latlong) {
        printf("could not init latlong\n");
        exit(1);
    }
    if (!pj_tgt) {
        printf("could not init target projection: %s\n", prj->init_string);
        exit(1);
    }

    printf("\tdescribe '%s'\n", prj->name);
    for (i=0; i<(sizeof(TEST_POINTS) / sizeof(point_t)); i++) {
        test_point=TEST_POINTS+i;
        x=test_point->x * DEG_TO_RAD;
        y=test_point->y * DEG_TO_RAD;

        pj_transform(pj_latlong, pj_tgt, 1, 1, &x, &y, 0);

        printf("\t\tit 'should forward transform (%f, %f) to (%f,%f)'\n", test_point->x, test_point->y, x, y);
        printf("\t\t\tvar prj=new Projections.%s();\n", prj->name);
        printf("\t\t\txy=prj.forward(%.20f, %.20f)\n", test_point->x, test_point->y);
        printf("\t\t\txy[0].should.equal_approximately %.20f\n", x);
        printf("\t\t\txy[1].should.equal_approximately %.20f\n", y);
        printf("\t\tend\n\n");


        x1=x;
        y1=y;
        pj_transform(pj_tgt, pj_latlong, 1, 1, &x1, &y1, 0);
        x1*=RAD_TO_DEG;
        y1*=RAD_TO_DEG;

        printf("\t\tit 'should inverse transform (%f, %f) to (%f,%f)'\n", x, y, x1, y1);
        printf("\t\t\tvar prj=new Projections.%s();\n", prj->name);
        printf("\t\t\txy=prj.inverse(%.20f, %.20f)\n", x, y);
        printf("\t\t\txy[0].should.equal_approximately %.20f\n", x1);
        printf("\t\t\txy[1].should.equal_approximately %.20f\n", y1);
        printf("\t\tend\n\n");
    }

    printf("\tend\n\n");
}

void projection::init_proj4() const
{
#ifdef MAPNIK_USE_PROJ4
    if (!proj_)
    {
#if PJ_VERSION >= 480
        proj_ctx_ = pj_ctx_alloc();
        proj_ = pj_init_plus_ctx(proj_ctx_, params_.c_str());
        if (!proj_)
        {
            if (proj_ctx_) {
                pj_ctx_free(proj_ctx_);
                proj_ctx_ = 0;
            }
            throw proj_init_error(params_);
        }
#else
        #if defined(MAPNIK_THREADSAFE)
        mapnik::scoped_lock lock(mutex_);
        #endif
        proj_ = pj_init_plus(params_.c_str());
        if (!proj_) throw proj_init_error(params_);
#endif
        is_geographic_ = pj_is_latlong(proj_) ? true : false;
    }
#endif
}

int        S57_initPROJ()
// NOTE: corrected for PROJ 4.6.0 ("datum=WGS84")
{
    if (FALSE == _doInit)
        return FALSE;

    const char *pj_ver = pj_get_release();
    if (NULL != pj_ver)
        PRINTF("PROJ4 VERSION: %s\n", pj_ver);

    // setup source projection
    if (!(_pjsrc = pj_init_plus(_argssrc))){
        PRINTF("error init src PROJ4\n");
        S57_donePROJ();
        return FALSE;
    }

    // FIXME: will need resetting for different projection
    _doInit = FALSE;

    if (NULL == _attList)
        _attList = g_string_new("");

    return TRUE;
}

int        S57_setMercPrj(double lat, double lon)
{
    // From: http://trac.osgeo.org/proj/wiki/GenParms (and other link from that page)
    // Note: For merc, PROJ.4 does not support a latitude of natural origin other than the equator (lat_0=0).
    // Note: true scale using the +lat_ts parameter, which is the latitude at which the scale is 1.
    // Note: +lon_wrap=180.0 convert clamp [-180..180] to clamp [0..360]

    const char *templ = "+proj=merc +lat_ts=%.6f +lon_0=%.6f +ellps=WGS84 +datum=WGS84 +unit=m";

    if (NULL != _pjstr) {
        PRINTF("WARNING: Merc projection allready set\n");
        return FALSE;
    }

    _pjstr = g_strdup_printf(templ, lat, lon);
    PRINTF("DEBUG: lat:%f, lon:%f [%s]\n", lat, lon, _pjstr);

    if (NULL != _pjdst)
        pj_free(_pjdst);

    if (!(_pjdst = pj_init_plus(_pjstr))) {
        PRINTF("ERROR: init pjdst PROJ4 (lat:%f) [%s]\n", lat, pj_strerrno(pj_errno));
        g_assert(0);
        return FALSE;
    }

    return TRUE;
}

// ===========================================================================
// method definitions
// ===========================================================================
GeoConvHelper::GeoConvHelper(const std::string& proj, const Position& offset,
                             const Boundary& orig, const Boundary& conv, int shift, bool inverse):
    myProjString(proj),
#ifdef HAVE_PROJ
    myProjection(0),
    myInverseProjection(0),
    myGeoProjection(0),
#endif
    myOffset(offset),
    myGeoScale(pow(10, (double) - shift)),
    myProjectionMethod(NONE),
    myUseInverseProjection(inverse),
    myOrigBoundary(orig),
    myConvBoundary(conv) {
    if (proj == "!") {
        myProjectionMethod = NONE;
    } else if (proj == "-") {
        myProjectionMethod = SIMPLE;
    } else if (proj == "UTM") {
        myProjectionMethod = UTM;
    } else if (proj == "DHDN") {
        myProjectionMethod = DHDN;
    } else if (proj == "DHDN_UTM") {
        myProjectionMethod = DHDN_UTM;
#ifdef HAVE_PROJ
    } else {
        myProjectionMethod = PROJ;
        myProjection = pj_init_plus(proj.c_str());
        if (myProjection == 0) {
            // !!! check pj_errno
            throw ProcessError("Could not build projection!");
        }
#endif
    }
}

void P4Projection::setDefinition(QString s)
{
    _definition = s;
    if (_pj != 0)
        delete _pj;
    _pj = pj_init_plus(_definition.toLatin1());
}

void project_inv(int *n, double *x, double *y, double *xlon, double *ylat, char **projarg){

  /* call the _inverse_ projection specified by the string projarg,
  * returning longitude and lat in xlon and ylat vectors, given the
  * numbers in x and y vectors (all vectors of length n) */

  int i;

  projUV p;
  projPJ pj;
  
  if (!(pj = pj_init_plus(*projarg)))
    error(pj_strerrno(*pj_get_errno_ref()));
/*  Rprintf("%s\n", pj_get_def(pj, 0));*/

  for(i=0;i<*n;i++){
    if(ISNAN(x[i]) || ISNAN(y[i])){
      xlon[i]=x[i];
      ylat[i]=y[i];
    } else {
      p.u=x[i];
      p.v=y[i];
      p = pj_inv(p, pj);
      if (p.u == HUGE_VAL || ISNAN(p.u)) {
	    Rprintf("inverse projected point not finite\n");
      }
      xlon[i]=p.u * RAD_TO_DEG;
      ylat[i]=p.v * RAD_TO_DEG;
    }
  }

  pj_free(pj);
}