I need to attribute a geographical area (based on shapefile polygons) to 8 million GPS points. There are 405 polygons. I currently have the code below, and a quick extrapolation tells me it will take 84 hours to compute this.
I am new to GIS, but intuitively, I feel that there has to be a smarter way around this than, for each point, randomly testing each polygon with a point-in-polygon function until there is match. For example, simply putting the points and polygons in 2 groups (such as "East" and "West") simply based on their easternmost/westernmost longitude may divide by two the number of operations to perform by almost 2.
Are there indeed such much efficient algorithms? Are any of them easily implementable in python?
from shapely.geometry import shape, Point, Polygon import pandas as pd from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap plt.figure(figsize=(16,16)) m = Basemap(projection='merc',llcrnrlat=15,urcrnrlat=55,llcrnrlon=70,urcrnrlon=140,lat_ts=20,resolution='c') m.readshapefile(china_shapefile, 'prefectures', drawbounds=True) print len(m.prefectures) #get data china_shapefile="C:/Users/.../CHN_adm2_SIMPLIFIED/CHN_adm2" points_csv="C:/Users/.../sample_map_data.csv" points=pd.read_csv(points_csv) points=points.dropna() lat=points["Latitude"].tolist() lon=points["Longitude"].tolist() coordinates=zip(lat,lon) for coordinate in coordinates: for info,shape in zip(m.prefectures_info,m.prefectures): poly = Polygon(shape) #print info["NAME_1"],info["NAME_2"] xpt,ypt = m(coordinate[1],coordinate[0]) point1 = Point(xpt,ypt) if point1.within(poly): print coordinate,info["NAME_1"],info["NAME_2"]#,info["NAME_3"] break
أكثر...
I am new to GIS, but intuitively, I feel that there has to be a smarter way around this than, for each point, randomly testing each polygon with a point-in-polygon function until there is match. For example, simply putting the points and polygons in 2 groups (such as "East" and "West") simply based on their easternmost/westernmost longitude may divide by two the number of operations to perform by almost 2.
Are there indeed such much efficient algorithms? Are any of them easily implementable in python?
from shapely.geometry import shape, Point, Polygon import pandas as pd from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap plt.figure(figsize=(16,16)) m = Basemap(projection='merc',llcrnrlat=15,urcrnrlat=55,llcrnrlon=70,urcrnrlon=140,lat_ts=20,resolution='c') m.readshapefile(china_shapefile, 'prefectures', drawbounds=True) print len(m.prefectures) #get data china_shapefile="C:/Users/.../CHN_adm2_SIMPLIFIED/CHN_adm2" points_csv="C:/Users/.../sample_map_data.csv" points=pd.read_csv(points_csv) points=points.dropna() lat=points["Latitude"].tolist() lon=points["Longitude"].tolist() coordinates=zip(lat,lon) for coordinate in coordinates: for info,shape in zip(m.prefectures_info,m.prefectures): poly = Polygon(shape) #print info["NAME_1"],info["NAME_2"] xpt,ypt = m(coordinate[1],coordinate[0]) point1 = Point(xpt,ypt) if point1.within(poly): print coordinate,info["NAME_1"],info["NAME_2"]#,info["NAME_3"] break
أكثر...