CODEWITHSUNDEEP
pythonshapefile
emax
emax

Reputation: 7255

Python: how to check if shape areas contain points?

I have a shapefile of Singapore that I visualize in this way:

import shapefile
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
### Plot shapefile
sf = shapefile.Reader("myShape.shp")
recs    = sf.records()
shapes  = sf.shapes()
Nshp    = len(shapes)
cns     = []
for nshp in xrange(Nshp):
    cns.append(recs[nshp][1])
cns = array(cns)
cm    = get_cmap('Dark2')
cccol = cm(1.*arange(Nshp)/Nshp)
#   -- plot --
fig,ax = plt.subplots(figsize=(20,10))

for nshp in xrange(Nshp):
    ptchs   = []
    pts     = array(shapes[nshp].points)
    prt     = shapes[nshp].parts
    par     = list(prt) + [pts.shape[0]]
    for pij in xrange(len(prt)):
        ptchs.append(Polygon(pts[par[pij]:par[pij+1]]))
        ax.add_collection(PatchCollection(ptchs,facecolor=cccol[nshp,:],edgecolor='k', linewidths=.1))     
ax.set_xlim(103.6,104.1)
ax.set_ylim(1.15,1.48)

enter image description here

Now I have a list of points with coordinates in a dataframe and I want to check in which region of the shapefile the points are.

mypoints
    lon         lat
0   103.619740  1.280485
1   103.622632  1.268944
2   103.622632  1.274714
3   103.622632  1.277600
4   103.622632  1.280485

In particular from the shapefile I can extract information of the specific area such as the name. In fact, in recs I have this information. For instance:

recs[0]:

[42,
 1,
 'STRAITS VIEW',
 'SVSZ01',
 'Y',
 'STRAITS VIEW',
 'SV',
 'CENTRAL REGION',
 'CR',
 '21451799CA1AB6EF',
 datetime.date(2016, 5, 11),
 30832.9017,
 28194.0843,
 5277.76082729,
 1127297.23737]

In this case STRAITS VIEW is the name of the area. At the the end I would like to have a dataframe like:

mypoints
    lon         lat       name
0   103.619740  1.280485  STRAITS VIEW
1   103.622632  1.268944  STRAITS VIEW
2   103.622632  1.274714  DOWNTOWN
3   103.622632  1.277600  BEDOK
4   103.622632  1.280485  CHANGI

Upvotes: 1

Views: 2222

Answers (1)

DarkCygnus
DarkCygnus

Reputation: 7838

To check if a coordinate lies inside some Polygon you could use the Intersects method from GDAL. This method compares two OGRGeometry objects and return a boolean value indicating whether they intercept or not.

With this you can iterate over your points, and for each one test if it Intersects() with all your Polygon areas. It should look something like:

for point_geom in your_points:
    #iterate over your area polygons and check 
    for area_geom in your_areas:
        if area_geom.Intersects(point_geom):
            #then you have a match, save or do as you want

Note: In case your points are not an OGRGeometry (your Polygons most probably are as you are reading them from a Shapefile), you can create such Geometry as suggested here:

from osgeo import ogr
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(1198054.34, 648493.09)

and in the case of Polygons:

ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(1179091.1646903288, 712782.8838459781)
ring.AddPoint(1161053.0218226474, 667456.2684348812)
#etc...add all your points, in a loop would be better
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

Upvotes: 1

Related Questions