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),float(sys.argv)) #check(-95,47)