5

I have a ESRI shapefile (from here: http://pubs.usgs.gov/ds/425/). I looking to use python to lookup information from the shape file (surficial material in this case) at a given latitude/longitude.

What is the best way to go about solving this problem?

Thanks.

Final solution:

#!/usr/bin/python

from osgeo import ogr, osr

dataset = ogr.Open('./USGS_DS_425_SHAPES/Surficial_materials.shp')
layer = dataset.GetLayerByIndex(0)
layer.ResetReading()

# Location for New Orleans: 29.98 N, -90.25 E
point = ogr.CreateGeometryFromWkt("POINT(-90.25 29.98)")

# Transform the point into the specified coordinate system from WGS84
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)
coordTransform = osr.CoordinateTransformation(
        spatialRef, layer.GetSpatialRef())

point.Transform(coordTransform)

for feature in layer:
    if feature.GetGeometryRef().Contains(point):
        break

for i in range(feature.GetFieldCount()):
    print feature.GetField(i)
arkottke
  • 547
  • 6
  • 14
  • If the point is not found in any of the `features`, then the last features will be treated as the match (and its fields printed). I would declare a separate variable `matched_feature` and assign to it immediately before the `break`, then use it for the next loop instead of the `feature` variable – Simon Ndunda Sep 25 '17 at 13:25

3 Answers3

3

Checkout the Python Shapefile Library

This should give you geometry and different info.

milkypostman
  • 2,955
  • 27
  • 23
  • Yes, it provides the geometry, but it doesn't provide the ability to see if a point is within a specified geometry. – arkottke Jan 04 '11 at 22:05
2

Another option is to use Shapely (a Python library based on GEOS, the engine for PostGIS) and Fiona (which is basically for reading/writing files):

import fiona
import shapely

with fiona.open("path/to/shapefile.shp") as fiona_collection:

    # In this case, we'll assume the shapefile only has one record/layer (e.g., the shapefile
    # is just for the borders of a single country, etc.).
    shapefile_record = fiona_collection.next()

    # Use Shapely to create the polygon
    shape = shapely.geometry.asShape( shapefile_record['geometry'] )

    point = shapely.geometry.Point(32.398516, -39.754028) # longitude, latitude

    # Alternative: if point.within(shape)
    if shape.contains(point):
        print "Found shape for point."

Note that doing point-in-polygon tests can be expensive if the polygon is large/complicated (e.g., shapefiles for some countries with extremely irregular coastlines). In some cases it can help to use bounding boxes to quickly rule things out before doing the more intensive test:

minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)

if bounding_box.contains(point):
    ...

Lastly, keep in mind that it takes some time to load and parse large/irregular shapefiles (unfortunately, those types of polygons are often expensive to keep in memory, too).

clint
  • 14,402
  • 12
  • 70
  • 79
2

You can use the python bindings to the gdal/ogr toolkit. Here's an example:

from osgeo import ogr

ds = ogr.Open("somelayer.shp")
lyr = ds.GetLayerByName("somelayer")
lyr.ResetReading()

point = ogr.CreateGeometryFromWkt("POINT(4 5)")

for feat in lyr:
    geom = feat.GetGeometryRef()
    if geom.Contains(point):
        sm = feat.GetField(feat.GetFieldIndex("surface_material"))
        # do stuff...
albertov
  • 2,314
  • 20
  • 15
  • I am having some difficulty specifying the point because there needs to be a coordinate transform from latitude/longitude to the reference system used by the shapefile. I assume this is somewhere in the gdal toolkit, but I can't seem to find it... yet. – arkottke Jan 04 '11 at 23:04
  • Look into osgeo.osr.SpatialReference and geom.Transform() (IIRC). I will update the example tomorrow (afk) – albertov Jan 04 '11 at 23:25
  • Yeah, I found some help on IRC and guided me through the transform. I updated the OP to reflect this information. Thanks for your help. – arkottke Jan 05 '11 at 00:26