I'm having trouble getting accurate mile markers through route traversal using PostGIS.
Based on the discussion in Geocode from mile marker/milepost to latitude/longitude?, I downloaded Texas Department of Transportation's official ShapeFile (link, see TxDOT Roadways 2010 at right).
I first imported the ShapeFile into Postgres using shp2pgsql. Then, using a C# program, I used PostGIS's ST_Line_Interpolate_Point function as the meat of my code that traverses each line one mile at a time. Here's the main part of the code:
public SortedDictionary<int, string> GenerateMileMarkers()
{
// make the sorted dictionary
SortedDictionary<int, string> dict = new SortedDictionary<int,string>();
// get length of route
double length = (double)(new NpgsqlCommand(
"SELECT ST_Length('" + geometry + "')",
connection)).ExecuteScalar();
// find each mile marker
// note that "increment" is 1/69.17, the reciprocal of miles/degree
for (int i = 0; i <= length/increment; i++)
{
double segmentLength = (i * increment)/length;
string mileMarkerPosition = (string)(new NpgsqlCommand(
// have to use ST_Linemerge to get from MULTILINESTRING to LINESTRING
// since ST_Line_Interpolate_Point requires LINESTRING
"SELECT ST_AsEWKT(ST_Line_Interpolate_Point(ST_Linemerge('" + geometry + "'), " + segmentLength + ")) As foo;",
connection)).ExecuteScalar();
dict.Add(i, mileMarkerPosition);
}
return dict;
}
I believe the SHP's reference measurement is a degree of latitude, so I did come calculations and figured that a degree has 69.17 miles, although thanks to Nicklas Avén's comment below, I now realize this is not correct since longitudinal lines are closer together the further north one goes.
Long story short, I generated serially-numbered points for each route that represent mile markers. The problem is I am coming up with significant errors.
I found the error by comparing Google Maps-indicated mile markers with those at a sampling of Interstate highway (IH) intersections or IH intersections with state borders. I'm finding error between ~4% and ~22. Most errors are approx ~15% +/- 1%. Here's an example of the error: where I-30 hits the Texas/Arkansas border is officially mile marker 223 (per Google Maps), but through this method I got 260. That's an error of about 16.6%. Also, where I-30 hits I-635 in Mesquite is about mile marker 56, but I calculated 67. That's a 19.6% error.
This variable error is really concerning to me. Since the errors are inconsistent, even within the same highway, it's not like I can apply a scaling factor to everything.
This is a problem because I am trying to do analysis on a few million events, and these events are expressed as happening at milepost X on highway Y.
Below a user rightfully questioned the projection. It appears the file may in fact have a projection, as it has a PRJ file:
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]]
Is it possible to make accurate measurements from this TxDOT-issued SHP?
(NOTE: I used "linear referencing" above, but I don't think this is true linear referencing because I am not correlating a point to a location along a route. The ST_Line_Interpolate_Point PostGIS function just happens to be in the Linear Referencing section of the PostGIS Reference.)
