Coordinate System C# – How to Convert Local Coordinates to WGS84 in C#

ccoordinate systemproj

Some people working on a the construction site use a local custom coordinate system. I have been given a proj4 string like this:

+proj=omerc +lat_0=46.325454996 +lonc=7.99050231 +alpha=60.4026421358 +gamma=90 +k=1 +x_0=350 +y_0=200 +datum=WGS84 +units=m +no_defs +type=crs

And x,y coordinates like this 429.079425606807, 261.126436507318

And I need to convert their coordinates into WGS84 and back

I tried to use the DotSpatial to make the transformation. Any other libraries in C# can be accepted.

 string proj4_custom = "+proj=omerc +lat_0=46.325454996 +lonc=7.99050231 +alpha=60.4026421358 +gamma=90 +k=1 +x_0=350 +y_0=200 +datum=WGS84 +units=m +no_defs +type=crs";

 DotSpatial.Projections.ProjectionInfo src = DotSpatial.Projections.ProjectionInfo.FromProj4String(proj4_custom);
 int googleEPSGCode = 4326;
 DotSpatial.Projections.ProjectionInfo trg = DotSpatial.Projections.ProjectionInfo.FromEpsgCode(googleEPSGCode);

 double[] xy = new[] { 429.079425606807, 261.126436507318 };
 double[] z = new[] { 0d };
 DotSpatial.Projections.Reproject.ReprojectPoints(xy, z, src, trg, 0, 1);

But the result in xy is not correct, on the map it's about 50 to 100 meters away. I should have a result really close to this correct result

Local coordinates: 429.07942560 261.1264365
WGS84: 7.99100318 46.326284519

What am I doing wrong?? If the problem is related to a bug in the libraries, what other libraries can I use that do it correctly?

Best Answer

I finally used the nuget package MaxRev.Gdal.WindowsRuntime.Minimal.

string proj = "+proj=omerc +lat_0=46.325454996 +lonc=7.99050231+alpha=60.4026421358 +gamma=90 +k=1 +x_0=350 +y_0=200 +datum=WGS84 +units=m +no_defs +type=crs";  SpatialReference local = new SpatialReference("");
local.ImportFromProj4(proj);
    
SpatialReference wgs84Reference = new SpatialReference("");
wgs84Reference.ImportFromProj4("+proj=latlong +datum=WGS84 +no_defs");
    
CoordinateTransformation coordinateTransform = new CoordinateTransformation(wgs84Reference, local);
double[] xy = new[] { lon, lat };
double[] z = new[] { alt };
coordinateTransform.TransformPoint(xy);
lon = xy[0];
lat = xy[1];

Related Solutions

[GIS] Slight translation with coordinate-transformation from Gauss-Kruger to WGS 84. Possible Reasons

Your Gauss-Krueger projection uses +datum=potsdam. Up to 2012, this was hard coded in proj4 to a very unprecise value using a 3-parameter-transformation.

You find more exact values for 7-parameter transformations in this topic:

http://forum.openstreetmap.org/viewtopic.php?id=12723

There is an even better ntv2-grid transformation available here (take the binary), that has to be in the same folder as your application and data, unlsss you specify full pathnames.

To compare the different possible values, I made this test batch file:

echo epsg31467-epsg4326 >out.txt
cs2cs +init=epsg:31467 +to +init=epsg:4326 31467.txt >>out.txt

echo proj-Definition epsg >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs  +to +init=epsg:4326 31467.txt >>out.txt

echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:4326 31467.txt >>out.txt

echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:4326 31467.txt >>out.txt

echo epsg31467-epsg3785 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3785 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3785 31467.txt >>out.txt

echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:3785 31467.txt >>out.txt

echo epsg31467-epsg3857 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3857 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3857 31467.txt >>out.txt

echo epsg31467-epsg900913 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:900913 31467.txt >>out.txt
echo epsg31467-proj900913 >>out.txt
cs2cs +init=epsg:31467 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt
echo proj31467-proj900913 >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt

with any sample Gauss-Krüger lon-lat coordinate pair in 31467.txt

[GIS] How to convert from a local coordinate system to UTM (and vice versa) using mathematical expression

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:

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):

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

Best Answer

Related Solutions

[GIS] Slight translation with coordinate-transformation from Gauss-Kruger to WGS 84. Possible Reasons

[GIS] How to convert from a local coordinate system to UTM (and vice versa) using mathematical expression

Related Question