4

I am now working with coordinate transformations between different CRSs, e.g. latlong to UTM.

After comparisons, I found out that the coordinate capture plugin in QGIS is quite accurate for the results.

Is there a way I can use the methods QGIS uses to do transformation in my own project?

I have checked that QGIS is using the proj4 library to do the transformation. I have implemented the proj4 transformation, but got slightly different results.

pj_svy21 = pj_init_plus("+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs") 

pj_latlong = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

double test1 = 37621.665;
double test2 = 12056.915;

projResult = pj_transform(pj_svy21, pj_latlong, 1, 1, &test1, &test2, NULL);

My result (latlong):
1.1253115030392806 103.9197679941208889

QGIS results:
1.35650722603 103.690059679 (check on map more accurate)

Can anyone help me on this?

Need to include iostream and proj_api.h:

int main () 

{

projPJ pj_svy21, pj_latlong

pj_svy21 = pj_init_plus("+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs")

pj_latlong = pj_init_plus("+proj=latlong +ellps=WGS84 +datum=WGS84 +no_defs")

double test1 = 12056.915;

double test2 = 37621.665;

projResult = pj_transform(pj_svy21, pj_latlong, 1, 0, &test1, &test2, NULL );

std::cout << "latlon: " << test1 * RAD_TO_DEG << ' ' << test2 * RAD_TO_DEG << std::endl;

return 0;

}
Improve this question

1 Answer 1

Reset to default
1

You may have swapped the coordinates somewhere.

If I reproject your resulting lonlat coordinates

103.919768  1.125312
103.690059679 1.35650722603

I get with cs2cs using:

cs2cs -f "%%.6f" +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs  inlon.txt >out.txt


37621.665436    12056.969913 0.000000
12056.914916    37621.664783 0.000000

Note that QGIS usualy expects X (horizontally or E) - Y (vertically or N) orientation.

Improve this answer
6
  • Hi AndreJ Thanks for the reply. I have changed the order of the north and east it seems ok but the output of pj_transform(pj_svy21, pj_latlong, 1, 1, &test1, &test2, NULL), test1 corresponds longitude while text2 corresponds latitude. Should it not be the other way around? Thanks
    – Ethan
    Commented Oct 24, 2014 at 14:43
  • No, you requested +proj=longlat, that means longitude first. The name of your variable pj_latlong is misleading.
    – AndreJ
    Commented Oct 25, 2014 at 7:19
  • ohh..i see..many thanks AndreJ, you solved my problems. :) Thanks.
    – Ethan
    Commented Oct 25, 2014 at 10:17
  • Hi AndreJ sorry to ask again, I have changed the +proj=longlat to +proj=latlong but the result are still long is first. I have pasted the code above you may try it out. Thanks in advance.
    – Ethan
    Commented Oct 25, 2014 at 11:08
  • For some strange reasons, +proj=latlong behaves the same way as +proj=longlat. Try +axis=neu, see trac.osgeo.org/proj/wiki/GenParms#Axisorientation
    – AndreJ
    Commented Oct 25, 2014 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.