Reading coordinates from shapefile using shapely

Question

I have coordinates and I am trying to find the roads they are on in a shapefile.

My code is as follows:

import shapefile
from shapely.geometry import Point # Point class
from shapely.geometry import shape # shape() is a function to convert geo objects through the interface


shp = shapefile.Reader(r'NBRoads/geonb_nbrn-rrnb_road-route.shp')


point_to_check = (-68.41856,47.28858) #example coordinate
all_shapes = shp.shapes() 
all_records = shp.records() 
print(shape(all_shapes[0]))    
for i in range(len(all_shapes)):
    boundary = shape(all_shapes[i]) 
    if boundary.distance(Point(point_to_check)) < 1e-8: 
       name = all_records[i][2] 
       print("The point is in", name) 

The problem is the LINESTRINGs coordinate system is not normal lat/long. LINESTRING (2433182.259100001 7613044.969799999, 2433256.279800002 7613043.196899999, 2433348 7613041, 2433368 7613040, 2433388 7613040, 2433411 7613039, 2433455 7613038, 2433500 7613037, 2433542 7613036, 2433609 7613034, 2433618.4628 7613033.746199999, 2433618.502300002 7613033.745099999) In my research I believe it has something to do with converting them using the .prj file from the shapefile but I cannot seem to make the link.

Here is the .prj info:

PROJCS["NAD_1983_CSRS_New_Brunswick_Stereographic",GEOGCS["GCS_North_American_1983_CSRS",DATUM["D_North_American_1983_CSRS",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",2500000.0],PARAMETER["False_Northing",7500000.0],PARAMETER["Central_Meridian",-66.5],PARAMETER["Scale_Factor",0.999912],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]

Answer 1

You can reproject your points from EPSG:4326 to EPSG:2953 (NAD_1983_CSRS_New_Brunswick_Stereographic) using pyproj like so:

import shapefile
import pyproj
from shapely.geometry import (Point, shape) # shape() is a function to convert geo objects through the interface
from shapely.ops import transform

from_crs = pyproj.CRS('EPSG:4326')
to_crs = pyproj.CRS('EPSG:2953')
transformer = pyproj.Transformer.from_crs(from_crs, to_crs, always_xy=True).transform

shp = shapefile.Reader('NBRoads/geonb_nbrn-rrnb_road-route.shp')
min_distance = 3  # distance where point is considered "on" a road

points_to_check = ((-68.41856,47.28858), (-68.4162, 47.28940)) #example coordinates

all_shapes = shp.shapes()
all_records = shp.records() 

for i in range(len(all_shapes)):
    boundary = shape(all_shapes[i])
    for lonlat in points_to_check:
        point_to_check = transform(transformer, Point(lonlat))
        if boundary.distance(point_to_check) < min_distance:
           name = all_records[i][2]
           print("The point is in", name)

If you have fiona installed, I prefer that to shapefile as you can read more data types than just shapefiles and the syntax is a bit cleaner:

import fiona
import pyproj
from shapely.geometry import (Point, shape)  # shape() is a function to convert geo objects through the interface
from shapely.ops import transform

from_crs = pyproj.CRS('EPSG:4326')
to_crs = pyproj.CRS('EPSG:2953')
transformer = pyproj.Transformer.from_crs(from_crs, to_crs, always_xy=True).transform

shapefile = 'NBRoads/geonb_nbrn-rrnb_road-route.shp'
min_distance = 3  # distance (in metres) where point is considered "on" a road
points_to_check = ((-68.41856, 47.28858), (-68.4162, 47.28940))  # example coordinates

with fiona.open(shapefile) as shp:
    for rec in shp:
        for lonlat in points_to_check:
            point_to_check = transform(transformer, Point(lonlat))
            boundary = shape(rec['geometry'])
            if boundary.distance(point_to_check) < min_distance:
                name = rec['properties']['name']
                print("The point is in", name)

You can also use geopandas which I like, but it has quite a different syntax. This example buffers the points, by the min distance, then uses a geopandas.sjoin spatial join:

import geopandas as gpd
from shapely.geometry import Point

from_crs = 'EPSG:4326'
to_crs = 'EPSG:2953'

shapefile = 'NBRoads/geonb_nbrn-rrnb_road-route.shp'
min_distance = 3  # distance (in metres) where point is considered "on" a road
points_to_check = ((-68.41856, 47.28858), (-68.4162, 47.28940))  # example coordinates

# separate the coordinates
lons, lats = list(zip(*points_to_check)) 

# create the geometry
points = [Point(lon, lat) for lon, lat in points_to_check]

# build the GeoDataframe
d = {'lon': lons, 'lat': lats,  'geometry': points}
gdfp = gpd.GeoDataFrame(d, crs=from_crs).to_crs(to_crs)

# Buffer the points
gdfp['geometry'] = gdfp.buffer(min_distance)

# Join lines to buffered points
gdfl = gpd.read_file(shapefile)
gdfp['name'] = gpd.sjoin(gdfp, gdfl, how='left')['name']

print(gdfp[['lon', 'lat', 'name']])

And if you have lots of points in a csv, you use the pandas read_csv method instead of looping through a list of points:

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

from_crs = 'EPSG:4326'
to_crs = 'EPSG:2953'

in_csv = 'test.csv'
results = 'test_results.csv'
shapefile = 'testline.shp'
min_distance = 3  # distance (in metres) where point is considered "on" a road

dfp = pd.read_csv(in_csv)  # add delimeter=', ' if csv has a space after the comma

# Build the GeoDataframe
gdfp = gpd.GeoDataFrame(dfp, crs=from_crs, geometry=gpd.points_from_xy(dfp.lon, dfp.lat))

# Reproject
gdfp = gdfp.to_crs(to_crs)

# Reproject and Buffer the points
gdfp['geometry'] = gdfp.buffer(min_distance)

# Join lines to buffered points with spatial join
gdfl = gpd.read_file(shapefile)
# Spatial join
gdfp['name'] = gpd.sjoin(gdfp, gdfl, how='left')['name']

# write out to csv
gdfp.replace(np.nan, '', regex=True)[['lon', 'lat', 'name']].to_csv(results)