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 = (-6841856,4728858) #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]]