4

I have two coordinate systems- UTM Zone 23S and a local one with the following conversion from UTM to Local:

EASTING: SQRT( (UTM_EAST-296787.355)^2 +UTM_NORTH-8095123.271)^2) *COS(ATAN((UTM_NORTH-8095123.271)/(UTM_EAST-296787.355))-0.0109372521874868)+6000

NORTHING: SQRT( (UTM_EAST-296787.355)^2+(UTM_NORTH-8095123.271)^2) *SIN(ATAN((UTM_NORTH-8095123.271)/(UTM_EAST-296787.355))-0.0109372521874868)+8000

I've spent a lot of time trying to get this to convert in proj4 code with QGIS, and the closest I can get is this:

+proj=omerc +y_0=8766 +x_0=6027.776 +alpha=-0.0109372521874868 +lat_0=-17.21294 +lonc=-46.91099472326864 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

which gets me within a few meters. But it is not perfect. The shapes do not align perfectly. Theres some rotation that I cannot figure out.

My question is this: is it possible to write a proj4 script to convert from one to the other using the above mathematical expressions?

  • Are you sure the UTM should be WGS84? Could be Aratu, Chua, Corrego Alegre, SAD69 or WGS72 as well. All those have a datum shift to WGS84. – AndreJ May 6 '16 at 18:19
  • To be honest I don't know. When using other local shapefiles the projection is only WGS84/Zone 23S. But the .prj files for some images show this: – joswhite May 6 '16 at 20:16
  • PROJCS["UTM_Zone_23_Southern_Hemisphere",GEOGCS["GCS_Geographic Coordinate System",DATUM["D_CORREGO_ALEGRE",SPHEROID["International_1909_Hayford_Intl_1924",6378388,297.0000000284015]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]] – joswhite May 6 '16 at 20:17
  • So this is Corrego Allegre. You might better use +ellps=intl +towgs84=-205.57,168.77,-4.12,0,0,0,0. For omerc, alpha is in degrees, while ATAN mostly delivers radians values. – AndreJ May 7 '16 at 5:38
2

The UTM 23S grid aligns with a degree grid only on the center longitude -45°. Other X and Y lines get bended (not only rotated) against a degree grid:

enter image description here

The oblique mercator projection is rotated against true north in the center point. The mathematical formula rotates the local grid against the UTM grid, which does not align with true north in the center point. So the omerc projection can not give identical values as the formula.

To get a good precision, I took the given center point from the formulas as Corrego Allegre UTM 23S (EPSG:22523), and reprojected it to Corrego Allegre degrees (EPSG:4225) to get the center coordinates. With that, I created the following custom CRS:

+proj=omerc +y_0=8000 +x_0=6000 +alpha=-0.061 +lat_0=-17.219614491 +lonc=-46.91098576 +gamma=0 +ellps=intl +towgs84=-205.57,168.77,-4.12,0,0,0,0 +units=m +no_defs

+towgs84 is taken from Corrego Allegre, and y_0 and x_0 from the formulas. From the formulas, I calculated a couple of points with LibreOffice Calc, and imported them as UTM 23S. The value +alpha is chosen empirically (for the reasons given above) to minimize the offset from the omerc grid to the points. It is about 0.4 meters for (10000 8000):

enter image description here

Note that the formula is only valid NorthEast of the center point. You see two of my points appearing at the "wrong" place.

  • @ AndreJ - Thank you! This is really good. Like you say, its not perfect but its within a meter everywhere that I have data. I spent so much time on this, working out the trig and banging my head on my desk. In the end, as you pointed out in the comments earlier, the corrego allegre was what I was missing. Also- Thank you for the diagrams. It was completely escaping me that only the center longitude line matched 'north' with the UTM grid. (Elementary I know, geographic vs projected coordinate systems). – joswhite May 9 '16 at 16:10
  • Also- @ ArMoraer, thank you for editing my question. As a native English speaker my grammar should really be better. – joswhite May 9 '16 at 16:17
  • Just going through your answer again- thanks for pointing out the degrees/radians distinction as well. – joswhite May 9 '16 at 16:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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