Skip to main content
2 of 3
edited tags
PolyGeo
  • 65.4k
  • 29
  • 114
  • 345

Reading coordinates from shapefile using shapely

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]]
Pearl
  • 113
  • 1
  • 4