I am working with a shapefile which can be found here. I am doing some data visualisation and want to see if certain longitude and latitudes are within a certain polygon in the shapefile.

This is the code I am using to check if a latitude and longitude is in a polygon. This works as I do it for other data

import geopandas as gpd
import shapely
from shapely.geometry import Point
from shapely.geometry import shape 
from shapely.geometry import Polygon
from descartes import PolygonPatch

//load in shapemap
soa_shape_map_path = './input/SOA2011_Esri_Shapefile_0/SOA2011.shp
soa_shape_map = gpd.read_file(soa_shape_map_path)

Then run through a set of longitude and latitudes in a dataframe to get what area it is part of

def check_df(df, polys, areas):
    # df contains crime data
    # polys contains all the paths dataframes containing the polygons
    # create a Polygon object from each of the dataframes

    for index, row in df.iterrows():
        lat = row["Latitude"]
        lng = row["Longitude"]
        point = Point((lng, lat))

    for row_idx, poly in enumerate(polys):
        if point.within(poly):
            df.at[index, "Areas"] = areas[row_idx]

return df

This never works because my longitude and latitude for df is always between 53/55 as it is Northern Ireland. When I check what the x and y co-ordinates are for soa_shape_map_path I get this

enter image description here

As you can tell, 54 is not even close to this number.

For my other shapefile map set I get this picture enter image description here

There definitely is a way to get the values as when I open the program, add the shp file in a layer and add an open street map it gives this

enter image description here

I have watched a couple YouTube videos and posts already but can not get what I need



  • Thanks Luke. I have literally just picked up this program today. Any idea on where to start? Nov 15 '18 at 2:12

Your soa_shape_map data is in a projected coordinate reference system (CRS) that uses metres for horizontal coordinate units (# TM65_Irish_Grid). You either need to project it to the same geographic CRS as your lon, lats (I don't know what geographic CRS you're using) or project your lon, lats to the same projected CRS as your shapefile.

E.g. (assuming WGS84 as the geographic CRS)

import geopandas as gpd

soa_shape_map_path = "path/to/SA2011_Esri_Shapefile_0/SA2011.shp"
soa_shape_map = gpd.read_file(soa_shape_map_path)

soa_shape_map_geo = soa_shape_map.to_crs(epsg=4326)  # EPSG 4326 = WGS84 = https://epsg.io/4326
print(soa_shape_map.iloc[0].geometry.centroid.y, soa_shape_map_geo.iloc[0].geometry.centroid.y)

377176.1212192121 54.62955604723917

Also, to avoid your explicit nested loops, you could have a look at geopandas spatial joins, e.g:

import geopandas as gpd

points_path = r"points.shp"
points_shape_map = gpd.read_file(points_path)
points_shape_map.to_crs(epsg=4326, inplace=True)

soa_shape_map_path = r"path/to/SA2011_Esri_Shapefile_0/SA2011.shp"
soa_shape_map = gpd.read_file(soa_shape_map_path)
soa_shape_map.to_crs(epsg=4326, inplace=True)

points_polys = gpd.sjoin(points_shape_map, soa_shape_map)  # Inner join (the default), only contains intersecting points (with the attributes of intersecting poly)
# points_polys = gpd.sjoin(points_shape_map, soa_shape_map, how="inner")  # Explicitly request an inner join
points_polys = gpd.sjoin(points_shape_map, soa_shape_map, how="left")  # Outer left join, contains all points and intersecting poly attributes (leave poly attribute None or NaN where there was no poly to intersect the point)
  • Thanks Luke. Just tried it and works perfectly. Will do some more reading about this topic. Nov 15 '18 at 2:31
  • @MartinO'Donnell also have a look at geopandas spatial joins to intersect the points (lon/lat dataframe) with your polygons. This removes the need for the explicit looping through each record, creating a Point object and then looping through each polygon record. It will likely be faster as well (if the polygons have a spatial index).
    – user2856
    Nov 15 '18 at 2:43
  • No problem. Will have a look in the morning. Thanks for your help Nov 15 '18 at 3:06
  • @MartinO'Donnell brief demo of poly to point spatial join in edit.
    – user2856
    Nov 15 '18 at 4:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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