OpenStreetMap

cs2cs help

Posted by h4ck3rm1k3 on 2 January 2010 in English.

I am doing the coord transform like that :

echo 358376.5 7753537 | cs2cs -v +proj=utm +south +ellps=intl +zone=24K +units=m -f "%.7f"a
# ---- From Coordinate System ----
#Universal Transverse Mercator (UTM)
# Cyl, Sph
# zone= south
# +proj=utm +south +ellps=intl +zone=24K +units=m
# ---- To Coordinate System ----
#Lat/long (Geodetic alias)
#
# +proj=latlong +ellps=intl
-40.3564136a -20.3105983a 0.0000000a

and it works, but now I want to do this in c++

see the code in MText.cpp
http://github.com/h4ck3rm1k3/TwoNickels

void convertPoint(double x, double y, double & rx, double & ry) 

{
projPJ pj_merc, pj_latlong;

if (!(pj_merc = pj_init_plus("+proj=utm +south +ellps=intl +zone=24K +units=m")) )
exit(1);

if (!(pj_latlong = pj_init_plus("+proj=latlong +ellps=intl")) )
exit(1);

point_offset, double *x, double *y, double *z );
double ax[1], ay[1], az[1];
ax[0] = x * DEG_TO_RAD;
ay[0] = y * DEG_TO_RAD;
az[0] = 0;

/** end of "caution" section. */

pj_transform(pj_merc, pj_latlong, 1, 1, ax, ay, az);

//printf("%.4f\t%.4f -> %.4f\t%.4f\n", *lat, *lon, y[0], x[0]);
rx = ay[0];
ry = ax[0];

}

Program produces ERROR::
node id='-905' version='1' ' lat='-26.92552447' ' lon='0.01349053807'
k='easting' v=358360.6875
k='northing' v=7753529

Please help!

mike

Discussion

Comment from greencaps on 3 January 2010 at 21:01

You make it impossible to help.

First the file MText.cpp is nowhere under the TwoNickels link.
Then you give the code for function convertPoint() but you do not show the parameter values when you call it.
You are not telling on which statement the error occurs.
And what is the error message? What has lat lon to do with calling convertPoint() ?

Comment from h4ck3rm1k3 on 4 January 2010 at 07:40

Sorry for the confusion.

MText is the new part of the Dime (two nickels)
http://github.com/h4ck3rm1k3/TwoNickels/blob/master/src/entities/MText.cpp

I have made a simple proj tool in dime as well here, but I forget to add that to git.
http://github.com/h4ck3rm1k3/TwoNickels/tree/master/proj/

On of my problems that was in the old code was the radians, I had converted the northing and easting to radians before converting. Now I do that after I get the results.

In fact, It works with the standard proj lib now, after I made the needed changes. Here is the current code :
So I dont *need* the proj patches, all the better.

Still, the code changes to proj are good, there were alot of bad typecasts, removing the constant from chars etc.

Ideally we could have a set of templates, one for each projection and do the code inline. It would be much faster.

to answer your questions : I had posted the parameters and results here with the code :

The new code using the standard proj interface :

void convertPoint(double x, double y, double & rx, double & ry) 

{
projPJ fromProj, toProj;
if (!(fromProj = pj_init_plus("+proj=utm +south +ellps=intl +zone=24K +units=m -f \"%.7f\" ")) )
exit(1);
if (!(toProj = pj_latlong_from_proj( fromProj )))
exit(1);
int islatlon= pj_is_latlong( fromProj );
int isgeocent=pj_is_geocent( fromProj );
char * name = pj_get_def( fromProj , 0);
std::cerr << "FROm NAME:"<< name << "|" << islatlon << "|" << isgeocent << std::endl;

islatlon= pj_is_latlong( toProj );
isgeocent=pj_is_geocent( toProj );
name = pj_get_def( toProj , 0);
std::cerr << "TO NAME:"<< name << "|" << islatlon << "|" << isgeocent << std::endl;

std::cerr.precision(10);
std::cerr.width(10);
std::cerr << x << "|" << y << std::endl;
// const double D2R =DEG_TO_RAD;
double ax = x;// * DEG_TO_RAD;
double ay = y;// * DEG_TO_RAD ;
double az = 0;

/** end of "caution" section. */
std::cerr.precision(10);
std::cerr.width(10);
std::cerr << ax << "|" << ay << "|" << std::endl;
int ret = pj_transform(fromProj, toProj,1,1, &ax, &ay, &az);

std::cerr.precision(10);
std::cerr.width(10);
std::cerr << ax << "|" << ay << "|" << ret << "|" << std::endl;

ax *= RAD_TO_DEG;
ay *= RAD_TO_DEG;

std::cerr << ax << "|" << ay << "|" << ret << "|" << std::endl;
rx = ay;
ry = ax;

pj_free(fromProj);
pj_free(toProj);

}


The inputs :
358376.5 7753537

the outputs :
FROm NAME: +proj=utm +south +ellps=intl +zone=24K +units=m|0|0
TO NAME: +proj=latlong +ellps=intl|1|0
358376.5|7753537
358376.5|7753537|
-0.7043522911|-0.3544868135|0|
-40.35641357|-20.31059831|0|
-20.3106:-40.3564

the convert point goes from a UTM northing/easting to lat/lon.
mike

Comment from greencaps on 4 January 2010 at 12:29

Sorry but I still see no parameters for

void convertPoint(double x, double y, double & rx, double & ry)

How did you call it was my question. Like:

convertPoint ( x, y, rx, ry );

Comment from greencaps on 4 January 2010 at 12:33

Sorry but I still see no parameters for

void convertPoint(double x, double y, double & rx, double & ry)

How did you call it was my question. Like:

double x = 12.23;
double y = 34.56;
double rx;
double ry;

convertPoint ( x, y, rx, ry );

and the result would be values in rx and ry

and the error would occur in that function.

You are showing all kinds of in and output that for me
has no connection to that functin.

Comment from h4ck3rm1k3 on 4 January 2010 at 12:48

Ahh,
Ok. Well the input and output are just

x=358376.5
y=7753537

outputs are :
ry=-20.3106
rx=-40.3564

see the main routine...
http://github.com/h4ck3rm1k3/TwoNickels/blob/master/proj/proj.cpp

it is all working now, so dont worry about it.

Log in to leave a comment