Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
add comment

2 Answers

up vote 2 down vote accepted

without looking thrugh all details in your problem, ther is one important misstake you are doing. You say that one degree is 69.17 miles. that is not right. one degree in east-west direction is different on different places. At the equator one degree is 1/360 * length of equator. at the poles one degree is zero lenght. But the longetudes are constant. I guess that is why you think you get incinsistent errors.

if you are not working over a too big area you should transform your data to some plannaer projection using miles or meters.

HTH

Nicklas

share|improve this answer
    
Thanks. I fixed everything by switching to 3081, which projects to meters. This was a good lesson for me! –  Aren Cambre Dec 28 '10 at 16:21
    
Yes, I saw you had great answers at the PostGIS list :-) –  Nicklas Avén Dec 28 '10 at 17:07
    
To follow that guidance peaople can look at: postgis.refractions.net/pipermail/postgis-users/2010-December/… , and press "next message" in the bottom to get all of it –  Nicklas Avén Dec 28 '10 at 17:09
add comment

Not only does the lat lon distance change with each degree, but it is not designed as a geodetic measurement system.

You should use a local coordinate system even if it is projected on the fly.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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