CODEWITHSUNDEEP
pythonshapely
Jonathan Lam
Jonathan Lam

Reputation: 1330

How to know if a coordinate is within a polygon in shape file not working?

I have coordinates (lat and lng) in an excel document. I have a shapefile which contains all the different Canada provinces shapes. I would like to be able to generate a new field in excel in order to classify the different coordinates into the different Canada provinces. I tried the below code, but it is not working.

import fiona
import shapely.geometry 

with fiona.open(r"D:\Users\Jonathan\Desktop\CRA-Project v2\Census Division\lcd_000b16a_e.shp") as fiona_collection:

shapefile_record = fiona_collection.next()

# Use Shapely to create the polygon
shape = shapely.geometry.asShape(shapefile_record['geometry'])
#print(shape)

point = shapely.geometry.Point(46.362914,-63.503809) # longitude, latitude

# Alternative: if point.within(shape)
if shape.contains(point):
    print("Found shape for point.")

Update 1: point = shapely.geometry.Point(46.362914,-63.503809)

Polygon: Link

Update Solution: I am updating this post because I found the solution and hopefully it can help someone!

  1. If you are using shapefiles, there is a projection associated with it (documentation)
  2. Review the documentation and identify which projection it is using. For example: In Canada it is Lambert conformal conic.
  3. Convert your lat/long to the projection equivalency
  4. Loop through the shapefile to identify if the new lat/long equivalent is within a polygon/multipolygon. You can do this with python!

Upvotes: 0

Views: 1452

Answers (1)

ewcz
ewcz

Reputation: 13087

There are several issues with this. First, your shapefile (Multipolygon) seems to be in a different projection so in order to work in longitude/latitude coordinates, one needs to transform it. Moreover, even after the transformation, the x-coordinate is the longitude, thus the coordinates of the point you are testing against should be swapped. An example might look as shown below. The input projection is assumed to be PCS_Lambert_Conformal_Conic (EPSG:3347) - this is specified in the accompanying prj file.

from functools import partial
import sys

import fiona
from shapely.geometry import Point, Polygon, asShape
from shapely.ops import transform
from shapely.wkt import loads
import pyproj

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:3347'),
    pyproj.Proj(init='epsg:4326'))

P = Point(-63.503809, 46.362914)
with fiona.open(sys.argv[1]) as F:
    for idx,feature in enumerate(F):
        G = transform(project, asShape(feature['geometry']))
        if G.contains(P):
            print(feature['properties'])
            break

This produces:

OrderedDict([('CDUID', '1102'), ('CDNAME', 'Queens'), ('CDTYPE', 'CTY'), ('PRUID', '11'), ('PRNAME', 'Prince Edward Island / Île-du-Prince-Édouard')])

i.e., it finds the coordinates within the Prince Edward Island.

Upvotes: 1

Related Questions