Possible issue with coordinate transform

Topics: General Topics
Oct 4, 2006 at 6:09 AM
I'm converting from the coordinate system given in my shapefiles to Lat/Long, and I'm coming up with obviously wrong answers. I don't know if it is the way the PRJ is being read by SharpMap, or if the manner in which I am transforming them is wrong.

ICoordinateTransformation trans2 = new CoordinateTransformationFactory()
.CreateFromCoordinateSystems(roads_shp.CoordinateSystem, GeographicCoordinateSystem.WGS84);
Point p =
trans2.MathTransform.Transform(new Point(2346961, 157068));

// Expected result of p is (-77.852023142309321, 34.1765389513421) (GPS confirmed)
// Actual result of p is (-75.1929603900457, 35.1059381758052)

The road_shp has a PRJ that as far as I know is correct:

PROJCS["NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
SPHEROID"GRS_1980",6378137.0,298.257222101
],
PRIMEM"Greenwich",0.0,
UNIT"Degree",0.0174532925199433
],
PROJECTION"Lambert_Conformal_Conic",
PARAMETER"False_Easting",2000000.002616666,
PARAMETER"False_Northing",0.0,
PARAMETER"Central_Meridian",-79.0,
PARAMETER"Standard_Parallel_1",34.33333333333334,
PARAMETER"Standard_Parallel_2",36.16666666666666,
PARAMETER"Latitude_Of_Origin",33.75,
UNIT"Foot_US",0.3048006096012192
]
Developer
Oct 4, 2006 at 7:18 AM
You are converting from NAD83 to WGS84 without specifying the TOWGS84 parameters, thus no datum-transformation is performed.
Oct 4, 2006 at 1:52 PM
Does that come from my PRJ and I'm missing it? Or does it go in the code? I haven't exactly seen any examples similar to this, and I'm not much of a GIS expert.
Oct 4, 2006 at 9:07 PM
I did the Lambert projection myself and I still do not know why the SharpMap code doesn't work as is. The calculated Eastings and Northings from the by hand projection are only off by 5 feet or so. Using the WGS84 latitudes and longitudes appeared to work properly.

LatFalse 33.7500000000 0.589048623
LongFalse -79.0000000000 -1.378810109
Lat1P 34.3333333333 0.59922971
Lat2P 36.1666666667 0.631227413
Efalse 2000000.0026166600
Nfalse 0.0000000000
a 20925646.3572663000
1/f 298.2572210100
e 0.081819191
e^2 0.00669438

m1 0.826650956
m2 0.808246494
t1 0.529982603
t2 0.509705998
tf 0.536504166
t' 0.531738035
n 0.577170255
F 2.06617083
rF 30183023.65
r' 30027971.45
theta 0.011554712

RealLat 34.1760120000 0.596483935
RealLon -77.8529620000 -1.358790519

CalcE 2346956.836
CalcN 157056.706
RealE 2346961
RealN 157068
Coordinator
Oct 11, 2006 at 9:42 PM
5 feet is pretty close...

That is about what you'll get using a WGS84 datum for a NAD83 projection.

Why do you need such accuracy? The coordinate conversion algorithm used in SharpMap (Bursa-Wolf) is only good to about 1-2 meters, and this is good enough for most purposes. A more rigorous conversion would be needed if you need more accuracy.
Oct 12, 2006 at 3:49 PM
I found out the real issue. SharpMap does not take the units into consideration. The data from the map is in US Survey Feet, yet the conversion is from Degrees to Meters and Meters to Degrees. I've added in the code to detect the need to convert the units since SharpMap fails to do this.

I may submit a patch, but I don't have the time at the moment.
Oct 12, 2006 at 9:31 PM
Wow. I just hit the exact same thing today. Can you save me some time by posting in the discussion where you made you changes and a code snippet?
Coordinator
Oct 12, 2006 at 11:30 PM
Units are part of the datum... if you are using WGS84, the units are meters. If you are using NAD, for example, the units are international feet. Are you specifying this in your projection datum?
Oct 13, 2006 at 4:05 AM
The projection file specifically states that it uses US survey feet and says how to convert it. I don't see where I'd do this any more "specifically" than in the file itself. There aren't any obvious API calls to make this conversion automatic.

I guess since there was somebody who got bit by the bug other than myself I'll post a patch when I clean it up.
Oct 13, 2006 at 4:11 AM
dlowther:

Right now I have the LambertConformalConic class check the units it is provided before doing any of the MetersToDegrees / DegreesToMeters calls. There is probably a cleaner way of doing this, so I'll look into it. As a quick fix you can just do these two changes (until a more permanent fix is available):

In the constructor make the change:
this.falseEasting = falseeasting.Value * LinearUnit.USSurveyFoot.MetersPerUnit;
this.falseNorthing = falsenorthing.Value * LinearUnit.USSurveyFoot.MetersPerUnit;

In DegreesToMeters:
return new SharpMap.Geometries.Point(
(rh1 * Math.Sin(theta) + this._falseEasting)/LinearUnit.USSurveyFoot.MetersPerUnit,
(rh - rh1 * Math.Cos(theta) + this._falseNorthing)/LinearUnit.USSurveyFoot.MetersPerUnit);

In MetersToDegrees:
double dX = p.X * LinearUnit.USSurveyFoot.MetersPerUnit - this._falseEasting;
double dY = rh - p.Y * LinearUnit.USSurveyFoot.MetersPerUnit
+ this._falseNorthing;

Of course this ONLY works if the units supplied are USSurveyFoot in your projection, replace with whatever unit is specified.

Also I noticed SharpMap ignores the conversion given for the degrees to radians specified in the GEOGCS block. Even though I'm sure the conversion in SharpMap is correct as-is, it is probably best to regard the units specified in the projection file as being exactly what should be used rather than those defined by SharpMap.
Developer
Oct 13, 2006 at 8:06 AM
After taking a closer look, I discovered that it is in fact your WKT that is wrong.
The parameters for your GRS80 datum is given in meters, although you state that the units are in feet. Try changing this value to feet as well, and the formulas should be correct.
Coordinator
Oct 13, 2006 at 9:05 AM
Ah, I see it now... missed it at the end of the projection. I thought you were saying the SharpMap projection was off by 5 feet or so... which would be impossible if there were the defect of improper unit conversion. But I see now it was your handiwork.

Good catch on the defect. I'll copy it in and unit test it and commit, unless Christian gets there first. Is the shapefile publicly available for me to use in testing?

BTW, you claim not to be a GIS expert and you are doing hand-crafted Lambert projections? ;)

Finally, on the matter of the degrees-to-radians conversion, which is specified in the PRJ file, there is only one right value - PI is a constant, as are the number of degrees in a circle. The option of reading it from the file is more consistent with the strategy of reading values from the file, but could then introduce errors if the file is corrupt or due to possible carelessness with floating-points. Also, if the PRJ is hand crafted or exported by a source which doesn't include it, SharpMap would have to use the defined value anyway. The important thing is the name of the measure. I don't really have a good answer for this, and can't find much in the specs and guidance about using it, but I'd tend to rely on the right answer vs. reading from the file. The file is meant to be human readable as well as machine readable, and the numeric values in it probably supports a human actor rather than a machine.
Coordinator
Oct 13, 2006 at 9:08 AM
SharpGIS -

But NAD83/HARN uses a GRS80 spheroid... and meausrements are in feet. Shouldn't Sharpmap use the datum's and not the spheroid's unit?
Oct 13, 2006 at 2:18 PM
The .prj file I used to create my coordinate system came straight from the ESRI directory. That doesn't AT ALL mean that it is right, only that it is one potential de-facto standard that - in my opinio - should be supprted if it is as simple to handle for as it appears to be.

Thanks for all the help.
Oct 13, 2006 at 4:06 PM
SharpGIS: I won't be able to change the PRJ in production because that is what is exported by the Arc tools. I have to take what is given to me (yay government shapefiles). I was under the impression that the UNIT at the end only applied to the PARAMETER's and not to the datum itself. The SPHEROID should have its own UNIT if it is not in meters (from what I've read).

codekaizen: I found a few papers describing the conversion and did it myself to find out why it wasn't working. The unit issues cropped up immediately when I tried to do it by hand ;D I'm slowing gaining momentum and learning all I need to about GIS. I'm working with my local city government on a surface street routing algorithm that takes into account traffic light timings and traffic flow data. I need to display it somehow and SharpMap/GIS won.
Coordinator
Oct 13, 2006 at 4:50 PM
dlowther -

ESRI changed the format of their PRJ files a while back to conform to the Well-Known Text standard from the OpenGIS Consortium (OGC). This is the actual standard that everyone uses.

Of course, ESRI is a major player in OGC, so it would seem they would have compelling reason to get it right. I'd trust their value. The thing with SharpMap is that not everyone using it has ESRI... ;)
Coordinator
Oct 13, 2006 at 5:05 PM
tunkeymicket -

Actually, SPHEROID doesn't have a unit in WKT. The GRS80 spheroid is always expressed in meters: 6378137.0 meters equitorial radius with an inverse flattening ratio of 298.257222101. Since the WKT specification appears not to allow for units to be specified, SharpMap has to know the definition of the spheroids in the various datums to get the units it uses. This is then converted to the Datum units during transformation.
Oct 13, 2006 at 6:58 PM
Gotcha. So should SharpMap make the conversion of the units in the projection before the calls to MetersToDegrees / DegreesToMeters? Because it really isn't MetersToDegrees using NAD83, its MapUnitsToDegrees and DegreesToMapUnits.
Coordinator
Oct 18, 2006 at 9:51 PM
This discussion has been copied to Work Item 4519. You may wish to continue further discussion there.