Reading the record values in several specific locations of a shapefile

Question

Hello everyone!

I have tried the following code:from osgeo import ogr

LandUse = ogr.Open("C:\Python26\ArcGIS10.0\landuse shpfile (polygon)\landuse.shp")
lyr = LandUse.GetLayerByName("landuse data_ypan")
lyr.ResetReading()
point = ogr.CreateGeometryFromWkt("POINT(503234.931458 3443753.43622)")
for feat in lyr:
    geom = feat.GetGeometryRef()
    if geom.Contains(point):
        sm = feat.GetField(feat.GetFieldIndex("Land_Group"))
        print sm

I am wondering if there is the possibility to import variables X, Y in the place of the number coordinates (red coloured text) in order to perform a for loop among several x, y coordinates.

Thank you very much.

irinivozinaki · Answer

Thank you very much for your help... it finally worked...

The only change that I have done is in:

points = [(x1,y1), (x2, y2),...,(xn, yn)]

I changed that to the following:

points=[(x,y)]

because of the fact that the code reads in every grid cell the (x, y) pair and as a result these pairs are always changing.

Thank you again!

ArkadiuszMatoszka · Answer

Not sure what you mean by saying "import variables X Y". You want to prepare a list or read them from file or something.
If list is ok for you, then code below should work (assume that coords are floats). If you want to read them from file the simply prepare list from text file before for loop:
points = []
f = open('C:\coords.txt', 'r')
for line in f.readlines():
    x, y = [coord.strip() for coord in line.split(';')] #suppose coords are delimited with semicolons
    points.append((float(x), float(y)))
f.close()

from osgeo import ogr

LandUse = ogr.Open("C:\Python26\ArcGIS10.0\landuse shpfile (polygon)\landuse.shp")
lyr = LandUse.GetLayerByName("landuse data_ypan")
lyr.ResetReading()
points = [(x1,y1), (x2, y2),...,(xn, yn)]
for coords in points:
    point = ogr.CreateGeometryFromWkt("POINT(%f %f)" % coords)
    for feat in lyr:
        geom = feat.GetGeometryRef()
        if geom.Contains(point):
            sm = feat.GetField(feat.GetFieldIndex("Land_Group"))
            print sm

Cheers
Arek