Python – Extracting Longitude and Latitude from Shapefile

geopandaslatitude longitudepythonshapefileshapely

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

//imports
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)
soa_geo_df=soa_shape_map['geometry']
soa_labels=soa_shape_map['SOA_LABEL']
reset_index(soa_geo_df)
reset_index(soa_labels)
soa_polys=soa_geo_df.values.T.tolist()

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]
            print("ADDED")            

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

https://www.youtube.com/watch?v=XI1qwzWPN1E

http://help.fulcrumapp.com/misc/guides-leveraging-external-tools/adding-latitude-and-longitude-data-to-shapefiles

Best Answer

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)
Related Question