[GIS] How to extract latitude/longitude values for specific points from a shapefile

arcmappostgisshapefile

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.

Best Answer

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

Related Solutions

[GIS] Indian Latitude Longitude values

This is a nice problem for OGR's python interface. This is my attempt, using world.shp. Here's my code snippet:

from osgeo import ogr # Import ogr
import random # Random number generators etc.

g = ogr.Open ("world.shp") # Open the shapefile with all the world's countries
# Loop over the different features until we find INDIA.
for feat in g.GetLayer(0):
    if feat.GetFieldAsString("NAME").find ("INDIA") >= 0:
        break # We found India, break out of the loop. 
              # feat now holds the feature of INDIA

geom = feat.GetGeometryRef() # Get the geometry of Indina
( lon_min, lon_max, lat_min, lat_max ) = geom.GetEnvelope() # Now get the envelope

n_points = 0 # A counter
india_lon = [] # A list where we'll store longitudes within India
india_lat = [] # A list where we'll store latitudes within India
while n_points < 1000: # Until we reach our required number of points...
    # Draw random longitudes and latitudes within the envelope
    r_lon = random.random()*(lon_max - lon_min ) + lon_min
    r_lat = random.random()*(lat_max - lat_min ) + lat_min
    # Create a point geometry
    pt = ogr.Geometry( ogr.wkbPoint )
    pt.SetPoint_2D( 0, r_lon, r_lat )
    # Is this geometry within india?
    if pt.Within ( geom ):
        # yes! Store the latitude&longitude and update counter
        india_lon.append ( r_lon )
        india_lat.append ( r_lat )
        n_points = n_points + 1
# Save to a CSV file
fp=open("india.txt", 'w' )
fp.write( "Longitude, Latitude\n" )# Header
for (lon, lat) in zip (india_lon, india_lat):
    fp.write ("%f,%f\n" % ( lon, lat ) )

fp.close()

And here's how the result looks like: One thousand random points over India

Not the most efficient way of dealing with this, but nice to learn :)

[GIS] Extracting Longitude and Latitude from shapefile

Your soa_shape_map data is in a projected coordinate reference system (CRS) that uses metres for horizontal coordinate units (# TM65_Irish_Grid). You either need to project it to the same geographic CRS as your lon, lats (I don't know what geographic CRS you're using) or project your lon, lats to the same projected CRS as your shapefile.

E.g. (assuming WGS84 as the geographic CRS)

import geopandas as gpd

soa_shape_map_path = "path/to/SA2011_Esri_Shapefile_0/SA2011.shp"
soa_shape_map = gpd.read_file(soa_shape_map_path)

soa_shape_map_geo = soa_shape_map.to_crs(epsg=4326)  # EPSG 4326 = WGS84 = https://epsg.io/4326
print(soa_shape_map.iloc[0].geometry.centroid.y, soa_shape_map_geo.iloc[0].geometry.centroid.y)

377176.1212192121 54.62955604723917

Also, to avoid your explicit nested loops, you could have a look at geopandas spatial joins, e.g:

import geopandas as gpd

points_path = r"points.shp"
points_shape_map = gpd.read_file(points_path)
points_shape_map.to_crs(epsg=4326, inplace=True)

soa_shape_map_path = r"path/to/SA2011_Esri_Shapefile_0/SA2011.shp"
soa_shape_map = gpd.read_file(soa_shape_map_path)
soa_shape_map.to_crs(epsg=4326, inplace=True)

points_polys = gpd.sjoin(points_shape_map, soa_shape_map)  # Inner join (the default), only contains intersecting points (with the attributes of intersecting poly)
# points_polys = gpd.sjoin(points_shape_map, soa_shape_map, how="inner")  # Explicitly request an inner join
points_polys = gpd.sjoin(points_shape_map, soa_shape_map, how="left")  # Outer left join, contains all points and intersecting poly attributes (leave poly attribute None or NaN where there was no poly to intersect the point)

Best Answer

Related Solutions

[GIS] Indian Latitude Longitude values

[GIS] Extracting Longitude and Latitude from shapefile

Related Question