1

I am creating some type of command line tool that will take a coordinate point(X and Y) and return an attribute from a polygon (shapefile) that surrounds the coordinate. I found this

To get attribute data of polygon layer by giving Lat/Lon

but it deals with a SQL statement and does not return any information.

I am using shapefile but I can convert it to another format.

I havent tried anything but I was hoping someone would give me a starting point.

1

You can execute the query based on spatial information, just as the link you provided says, but as you note this will return a new shapefile that intersects that coordinate (not the attributes). You can take an added step in Python to return desired attributes for the newly created selection shapefile.

Via command line:

ogr2ogr -sql "SELECT geometry FROM [base shapefile name without *.shp, eg input] WHERE ST_Intersects(ST_GeomFromText('POINT([X] [Y])'), [input].geometry)" -dialect SQLITE [output.shp] [input.shp]

Then, as a Python command:

from osgeo import ogr
shapefile = [output.shp]
driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(shapefile, 0)
layer = ds.GetLayer()

for feature in layer:
    print feature.GetField("AttributeField 1")
    print feature.GetField("AttributeField 2")

etc. until you list all the fields you want to query for the shapefile selected from your original point. If you want to avoid going into Python altogether, you can either put that section in its own .py file and call it via command line with "Python file.py", or else write in command line "Python -c "from osgeo import ogr; shapefile="... and so on.

Here's a great resource for getting features from vector data with OGR and Python, if you don't mind going in that direction: https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html

  • Why do I get an unexpected identifier error when using ogr2ogr -sql "SELECT * FROM polygons WHERE ST_Contains ST_GeomFromText('POINT(1305365 406287)')" output.shp polygons.shp – GreyHippo Oct 25 '17 at 19:42
  • Wow, I made two mistakes - see the edited answer above. Dialect flag (-dialect SQLITE) needs to be set to use the ST_ functions, and the OID column needs to be set to read and copy out individual features from the individual features (to find the intersection with the point). – AlecZ Oct 25 '17 at 20:02
  • I tried the ogr2ogr command but it selects all features of the input shapefile ogr2ogr -sql "SELECT * FROM lakes WHERE ST_Intersects(ST_GeomFromText('POINT(-220475 3482745)'), cat)" -dialect SQLITE output4.shp lakes.shp I am using a sample dataset from alaska lakes. – GreyHippo Oct 25 '17 at 20:27
  • Ah, I see where the issue is. Selecting based on OBJECTID like I thought is making the intersect statement return everything. Try the modified code above, also printed here: – AlecZ Oct 25 '17 at 21:56
  • ogr2ogr -sql "SELECT geometry FROM [base shapefile name without *.shp, eg input] WHERE ST_Intersects(ST_GeomFromText('POINT([X] [Y])'), [input].geometry)" -dialect SQLITE [output.shp] [input.shp] – AlecZ Oct 25 '17 at 21:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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