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 have obtained public crash data for the state of Minnesota. I would like to be able to programmatically extract the latitude/longitude of each crash using this data in combination with a series of shapefiles that cover all roads in the state. For each crash, I believe I have enough information to do this, but my knowledge of GIS is keeping me from getting things done.

I'll start with a simple example in the hopes that someone can point me in the right direction. You'll have to forgive me if my GIS terminology is not what it should be. The shapefile I'm using is here: http://www.dot.state.mn.us/maps/gisbase/datafiles/county/ramsey.zip

Accident 090710194, for example, should be at mile 1.228 along TIS_CODE 0462000015. When I use a tool like Quantum GIS to find all attributes with a TISCODE of 0462000015, there are 81 features to choose from. If I select the feature with a BEGM (beginning mile) and ENDM (ending mile) that encapsulates 1.228, I find one feature that goes from BEGM of 1.21 to ENDM of 1.3. This puts me "in the ballpark" for a precise location. Then, I can use Quantum GIS' "Coordinate Capture" plugin, to guess that it's somewhere around -93.20228,45.05062. I know these measurements are right, because a double check on Google Maps shows the same approximate location as Quantum GIS: http://maps.google.com/maps?q=45.05062,+-93.20228

Really, there are 2 problems with the above approach: 1) It's very slow; I need to do this for over 500K crashes, so I need to do something programmatically, and 2) It's not very precise.

I've been told by folks in the Minnesota Department of Public Safety that if I had ESRI ArcMap (the tool they use), I could "create a point event", which sounds like a precise measurement along the feature I selected, and from there I would see the precise location, and (I'm hoping) could probably get latitude/longitude. I don't have this software, but it sounds something like what is shown in the following video, where toward the beginning of this video they're creating a "route event". Since I'm not allowed to post more than 2 links in this post, in order to see what I mean you'll need to Google "arcgis point event" and select the first link, which is to webhelp dot esri dot com.

Everything boils down to 2 questions for me right now: 1) Is programmatically extracting this information possible, even with ESRI ArcMap, and 2) If so, what (preferably open source) tools can aid me? PostGIS? Geotools? Geoserver? Any advice at all would be appreciated.

share|improve this question
add comment

3 Answers

You can use PostGIS to achieve what you need. PostGIS has a function called ST_line_interpolate_point(geom, location) which can take a line string e.g your road, and a location which is a percentage from the start of the road e.g distance/road_length.

I would, as I did this at work last week for something like this, create a table with your Accidents (events) and a table with your roads (just upload the one you have from QGIS into PostGIS using SPIT) and then join them using a common attribute.

SELECT ST_line_interpolate_point(roads.the_geom, event.distance / ST_Length(roads.the_geom)), event.id
FROM event
JOIN roads on event.roadid = roads."TIS_Code"

The above syntax might not be 100% but should get you some of the way. You might have to tweak it to handle if the road is shorter then the event distance etc (this might work added to your join WHERE event.distance > roads.BEGM AND event.distance < roads.ENDM)

If you want to then get the corrdinates from the points you can use the range of ST_As... functions or you can wrap the above query in a another query.

SELECT ST_X(point),ST_Y(point), id
FROM (
    SELECT ST_line_interpolate_point(roads.the_geom, event.distance / ST_Length(roads.the_geom)) as point, 
           event.id as id
    FROM event
    JOIN roads on event.roadid = roads."TIS_Code")

Works like a charm.

EDIT:

This might work better:

SELECT ST_line_interpolate_point(roads.the_geom, event.distance / ST_Length(roads.the_geom)), event.id
FROM event
JOIN roads on event.roadid = roads."TIS_Code" AND event.distance > roads.BEGM AND event.distance < roads.END
share|improve this answer
    
Thank you, this is exactly what I was looking for! I'll give it a shot in the next couple of days. –  user5570 Jan 22 '12 at 15:20
    
@Nathan, Do you need a sql statement, in your answer, to order the roads by BEGM so ST_line_interpolate_point always starts at the correct end of the Road? Or did you do this when you "join them using a common attribute" ? –  klewis Jan 22 '12 at 18:17
    
@klewis Not really. ST_line_interpolate_point will always start at the start of the line as per the direction that it was digitized, so ordering based on BEGM will have no affect. –  Nathan W Jan 22 '12 at 20:18
    
I get st_x = 484008.838340942, st_y = 4988550.07593788, which is unfortunately not latitude and longitude. However, in QGIS, I see this coordinate system below latitude and longitude (by an image of a grid, rather than the image of a globe for latitude/longitude). I feel that I am very close, is there another function that can convert these values for me? –  user5570 Jan 28 '12 at 5:31
    
To answer my own question, the answer was to convert to SRID 4326. Here is the final query that led me to success: select st_x(point), st_y(point), accn from (select st_line_interpolate_point(st_transform(LineMerge(roads.the_geom), 4326), events.distance / st_length(roads.the_geom)) as point, events.accn FROM events JOIN roads on events.road_id = roads.tis_code AND events.distance > roads.begm AND events.distance < roads.endm) as point; –  user5570 Jan 28 '12 at 5:40
add comment

Yes, you can extract this information with ESRI ArcMap. You won't need any extra tools besides what ArcMap provides. What you'll be doing is linear referencing. You will need to create a route in ArcMap which should not be difficult with the data you have. Each accident number and corresponding mile number will be in an event table. You can plot these events on the route using the "make route event layer" tool. You can then create a feature class from this and add fields for lat/lon.

share|improve this answer
    
Well, at least I know it's possible now. Unfortunately, the pricing of ArcGis (of which ArcMap is a part, I believe) is way outside my budget for this project: esri.com/software/arcgis/arcgis-for-desktop/pricing.html –  user5570 Jan 21 '12 at 17:51
    
You can always get a 60 day trial here with all the features you'd need: esri.com/software/arcgis/arcgis-for-desktop/free-trial.html –  tbuckler Jan 21 '12 at 20:10
add comment

Here is your dataset in Fusion Tables: Table: http://www.google.com/fusiontables/DataSource?snapid=S3683355IEK Map: http://www.google.com/fusiontables/DataSource?snapid=S368337_PGn

This doesn't solve your issue about extracting a point from each line segment, but it is a handy, free, and easy way to visualize the data on a map.

share|improve this answer
    
Not really a answer. This would be better as a comment on the main question. –  Nathan W Jan 23 '12 at 0:02
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.