Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Python: assign point to roads shapefile

I have some data points from gps data (Lat,Lon) and the relative shapefile of the area I am considering. So:

import pandas as pd
import shapefile
sf = shapefile.Reader('roadshape')
df = pd.read_csv('gpsdata.csv')

Now df contains the data I am analyzing is something like this:

 ID          x          y
3447    11.400427   48.816806
3448    11.400759   48.816772
3449    11.401424   48.816684
3450    11.401758   48.816631
3451    11.402090   48.816566
3452    11.402422   48.816490

I want to assign to each the respective shapefile segment. I am trying to do the following. I am considering the bounding box and I want to try to see in which bounding box the data points the data bar.

dfs = pd.DataFrame()
shapes = sf.shapes()
X =list()
Y=list()
for i in range(0,len(shapes)):
      X.append([shapes[i].bbox[0],shapes[i].bbox[2]])
      Y.append([shapes[i].bbox[1],shapes[i].bbox[3]])
dfs['X'] = X
dfs['Y'] = Y

Now how can I check in which bbox my points are? Where dfs is something like this

dfs = 
                 X                             Y
0   [10.9467244189, 10.9704393002]  [48.2671975178, 48.2697440003]
1   [11.5138847999, 11.5143541004]  [48.6497096997, 48.6515363002]
2   [11.4618209998, 11.4620896004]  [48.9305448001, 48.9307776004]
3   [10.6196591004, 10.6207268996]  [48.8635958001, 48.8665684003]
4   [10.652098, 10.6559025999]  [48.8005320998, 48.8042877999]
5   [11.1863882997, 11.1884544004]  [48.3726685999, 48.3738253996]
6   [11.1580075998, 11.1593822] [48.3785226999, 48.3791247996]
7   [11.1077987, 11.1112508996] [48.3829125003, 48.3830440999]
8   [11.0842697004, 11.0886483996]  [48.3840543003, 48.3879626001]
9   [11.0910959001, 11.0926532003]  [48.3903297003, 48.3916850002]
10  [11.4766434001, 11.4822778002]  [49.0389071001, 49.0399456003]
11  [11.7037148998, 11.7073818] [48.6927748996, 48.6961230001]
12  [11.7767894997, 11.7770049998]  [48.6279809001, 48.6279908997]
like image 243
emax Avatar asked Oct 25 '15 23:10

emax


Video Answer


1 Answers

I would go for fiona,shapely and built-in csv!

Say i have two files:1) A line shapefile with several line records in it. 2) A csv file with three columns i.e. id, longitude(x), latitude(y)

CSV file content (lat, long is UTM)-

id,lat,long

0,207726.012448,2733349.10914
0,197599.591396,2730510.17345
0,203187.5176,2736670.5686
0,207301.877268,2730639.81898
0,200929.610894,2726377.9799
0,204604.214301,2737608.342
0,203780.386032,2734372.2709
0,203077.172106,2731166.44271
0,202477.371994,2728622.46292
0,202249.861606,2734889.33996
0,201794.840831,2732159.21531

Now below code generates bbox for each record of the line shapefile and checks what point of csv file is inside of that bbox after all it prints status!

import fiona
from shapely import geometry
import csv
csv_file = open(r"C:\data_test.csv",'rb')
reader = csv.reader(csv_file)
pnts = []
for i in reader:
    pnts.append ((i[1],i[2]))
bboxes = []   
fiona_collection = fiona.open(r"C:\Rd.shp")
for i in fiona_collection:
    bboxes.append(geometry.asShape(i['geometry']).bounds)


for i in pnts:
    for j in bboxes:
        shape = geometry.box(*map(float,j))
        pnt = geometry.Point(*map(float,i))
        if shape.contains(pnt):
            print "bbox {0} contains {1} ".format(j,i)

It prints-

bbox (197306.68136754428, 2729718.6493367185, 197823.7504301681, 2732035.118737273) contains ('197599.591396', '2730510.17345') 
bbox (202799.29018572776, 2728807.5534000006, 204604.21430105247, 2737608.342) contains ('203187.5176', '2736670.5686') 
bbox (202799.29018572776, 2728807.5534000006, 204604.21430105247, 2737608.342) contains ('203780.386032', '2734372.2709') 
bbox (202799.29018572776, 2728807.5534000006, 204604.21430105247, 2737608.342) contains ('203077.172106', '2731166.44271') 
bbox (201237.0747999996, 2727558.2885, 203336.7932000002, 2728807.5534000006) contains ('202477.371994', '2728622.46292') 

you can download test data at here.

like image 102
SIslam Avatar answered Sep 19 '22 09:09

SIslam