2

I have a csv file with two sets of lat/lon points - the first set of points is represented by columns V1(lon) and V2(lat), and the second set of points by V3(lon) and V4(lat) (the csv does not have headers).

> head(merge)
                 V1               V2                V3               V4   V5     V6
1 -122.030277484219 38.0534546969437 -122.263688968328 38.0180951886275    0 Sunday
2 -122.012178329029 38.5883736394547 -121.269556181423 37.9550435479333    1 Sunday
3 -121.817242008376 38.8419619506953 -121.305704278949 37.2665873119052   10 Sunday
4 -122.100160967342 38.7339037710083 -121.839444904779 37.8529730771058  100 Sunday
5 -122.424897962683 37.3465530990606 -122.193097527866 37.2373493348652 1000 Sunday
6 -122.262026449704 37.6901105931389 -121.913111829827 37.0947943964057  101 Sunday

In python, I need to check the location of each point - and create a separate column in my .csv file containing the name of the feature in which the point is located. I have a shapefile containing polygons - lets call it counties.shp - so basically I need to find which polygon each point falls within.

I understand that first I have to import my shapefile, transform it to WGS84, import the csv, create points from lat lons, then loop over the points to search for overlapping features, although I am not sure about the syntax in python. I have found the following resource but it needs to be modified. Any ideas where to start?

#!/usr/bin/python
import ogr
from IPython import embed
import sys

drv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape file
ds_in = drv.Open("counties.shp")    #Get the contents of the shape file
lyr_in = ds_in.GetLayer(0)    #Get the shape file's first layer

#Put the title of the field you are interested in here
idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("P_Loc_Nm")

#If the latitude/longitude we're going to use is not in the projection
#of the shapefile, then we will get erroneous results.
#The following assumes that the latitude longitude is in WGS84
#This is identified by the number "4326", as in "EPSG:4326"
#We will create a transformation between this and the shapefile's
#project, whatever it may be
geo_ref = lyr_in.GetSpatialRef()
point_ref=ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)
ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref)

def check(lon, lat):
    #Transform incoming longitude/latitude to the shapefile's projection
    [lon,lat,z]=ctran.TransformPoint(lon,lat)

    #Create a point
    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat)

    #Set up a spatial filter such that the only features we see when we
    #loop through "lyr_in" are those which overlap the point defined above
    lyr_in.SetSpatialFilter(pt)

    #Loop through the overlapped features and display the field of interest
    for feat_in in lyr_in:
        print lon, lat, feat_in.GetFieldAsString(idx_reg)

#Take command-line input and do all this
check(float(sys.argv[1]),float(sys.argv[2]))
#check(-95,47)
  • You could start by modifying the last active line to loop through the file, use standard python file operations and then call check() for each line.. apart from that it looks fairly good (but I'd get rid of the Iron Python parts as I don't have that installed). On my install the syntax of the OGR import is 'from osgeo import ogr', try it in your IDE and see which one works for you. – Michael Stimson May 11 '16 at 2:24

Your Answer

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

Browse other questions tagged or ask your own question.