1

I have large set of GPS coordinates of human mobility in different parts of China, but I am only interested in points in Beijing City. These points are in excel csv file that I read into pandas dataframe using python 3.

I want to filter all GPS points inside Beijing city and discard all other points in the dataframe. One approach that I have in my mind is to define a rectangular bounding box and note down the following coordinates (latitude, longitudes) of bounding box.

Bottom Left  (x1,y1) = 39.77750000, 116.17944444  
Top Left     (x1,y2) = 40.04722222, 116.17944444  
Bottom Right (x2,y1) = 39.77750000, 116.58888889  
Top Right    (x2,y2) = 40.04722222, 116.58888889  

Then for each point, compare its latitude and longitude values with all the vertices using multiple if else conditions, but I think their should be more optimal approaches to solve this in simple way. Is there any easier way or supporting library to solve this problem?

Example image is attached for clear understanding for the problem

enter image description here

The input dataset in given in the link https://drive.google.com/file/d/1JjvS7igTmrtLA4E5Rs5D6tsdAXqzpYqX/view

1 Answer 1

3

Use within:

Returns a Series of dtype('bool') with value True for each geometry that is within other.

from shapely import geometry
import geopandas as gpd
import pandas as pd

file = r'C:\folder\user141.csv'
crs = {'init': 'epsg:4326'}

#Create a geoseries holding the single polygon. Coordinates in counter-clockwise order
pointList = [(39.77750000, 116.17944444),(39.77750000, 116.58888889),(40.04722222, 116.58888889),(40.04722222, 116.17944444)]
poly = geometry.Polygon(pointList)
spoly = gpd.GeoSeries([poly],crs=crs)

#Create geodataframe of points
dfcsv = pd.read_csv(file)
geometry = [geometry.Point(xy) for xy in zip(dfcsv.lat, dfcsv.long)]
dfpoints = gpd.GeoDataFrame(dfcsv, crs=crs, geometry=geometry)

#Create a subset dataframe of points within the polygon
subset = dfpoints[dfpoints.within(spoly.geometry.iloc[0])]
print('Number of points within polygon: ', subset.shape[0])

Outputs:

Number of points within polygon:  58772
2
  • Works great! Two things:
    – Jheasly
    Feb 16, 2022 at 23:06
  • In geometry = [geometry.Point(xy) for xy in zip(dfcsv.lat, dfcsv.long)] I think you mean dfcsv.lng (dfcsv.long throws an AttributeError). And FWIW, I was able to run it w/o the crs parameter.
    – Jheasly
    Feb 16, 2022 at 23:13

Not the answer you're looking for? Browse other questions tagged or ask your own question.