1

I have a huge dataset of points and roads. The point file is of Residential houses or Comercial shops/ markets. And the Line network is the roads around them. Ideally, I want to take 1 point (input) and find out its distance from nearest 4 lines (roads). The roads shapefile also includes some attributes (including name and Width). The final output should include Distance(meters) from inputted point to the nearest 4 roads and the attributes of these 4 roads.

roads = gpd.read_file (r'fRoads.shp')
long = float(input("Enter longitude E value ))
lat = float(input("Enter latitude N Value ))
radius = float(input("Enter Radius value in km: "))
p1 = Point((long, lat))
df = pd.DataFrame({'a':[lat,long]})
po = gpd.GeoDataFrame(geometry = [p1], crs = roads.crs)
bufo = createbuf(po,radius)
cliped = createClip(bufo,roads)
fig, ax= plt.subplots(figsize = (15,7))
roads.plot(ax=ax, color = 'black')
bufo.plot(ax=ax, color = 'None', edgecolor = 'red')
po.plot(ax=ax, color = 'blue')

fig, ax= plt.subplots(figsize = (15,7))
bufo.plot (ax = ax, color ='None', edgecolor = 'Black')
cliped. plot (ax = ax, color = 'red')
po.plot (ax = ax, color = 'Black')
print (cliped.length)
dist  = po.distance(Point (roads))
print (dist)

I am using the Distance Function from shapely. but it gives the distance between 2 points not between a point and line.

The Image contains a point and few lines. I want the Distance in meters from point to all lines

I need something equivalent to the following code but without Arcpy and that should be Open Source.

    import arcpy
    from openpyxl import load_workbook
    
    arcpy.env.workspace = r'C:\Users\MJ\Desktop\populationFunc\roads func\froads2'
    arcpy.env.overwriteOutput = True
    
    long = float(input("Enter longitude N value e.g: 77.xxxxx: "))
    lat = float(input("Enter latitude E Value e.g: 28.xxxxx:  "))
    radius = float(input("Enter Radius value in km: "))
    spatial_reference = arcpy.SpatialReference(32642)
    p1 = arcpy.Point(long, lat)
    pnt_geometry = arcpy.PointGeometry(p1, spatial_reference)
    roads = 'fRoads.shp'
    point = r'C:\Users\MJ\Desktop\populationFunc\roads func\propData\propPoints.shp'
    search_radius = '1500 Meters'
    location = 'NO_LOCATION'
    angle = 'NO_ANGLE'
    out_table = "out_table.dbf"
    closest = 'ALL'
    closest_count = 4
    arcpy.GenerateNearTable_analysis(pnt_geometry, roads, out_table, search_radius, location, angle, closest, closest_count)
    fp = arcpy.SearchCursor(out_table)
    for f in fp:
        print f.getValue("NEAR_FID")
        print f.getValue("NEAR_DIST")
4
  • Did you try to get the Points of the lines (lines consists of Points), calculate the distance of point to each one of these Points and keep the minimum one? Mar 22 at 6:27
  • 1
    Distance Function from shapely. gives also the distance between a point and a line
    – gene
    Mar 22 at 9:17
  • And the shapely distance is euclidean, so you can use it with longitude, latitude (angles)
    – gene
    Mar 22 at 9:33
  • Can you share any example... I have tried but could not get done with that. I have seen point-to-point distance but not point to line.
    – alauddin
    Mar 22 at 21:13
1

I think your issue is that your road geometries are using a projected coordinate system (UTM Zone 42N (EPSG Code 32642); according to your second code snippet with arcpy) with metre for the unit, while you are using lat/lon coordinates (in degrees) for your point. However, I was not able to run your script to double check.

If you need to use lat/lon coordinate input but have your data in a projected coordinate system, then you may be able to use shapely.ops.transform(func, geom). I have added an example at the end of this answer.

Since I was not able to get your code working as well as I have no data available in the coordinate system your were using, I rewrote what I think you are looking for, using shapely and geopandas. Note that I use another dataset and projection (UTM Zone 20N) than you, neglected the user input and use a hard-coded point coordinate and buffer radius.

import matplotlib.pyplot as plt
import geopandas as gpd
from shapely import geometry
from shapely.geometry import Point

# Using test data from here: https://nsgi.novascotia.ca/gdd/
# NAD83(CSRS) / UTM Zone 20N (EPSG:2961)
test_data = r"data\TRNS_NSRN_Addressed_Roads_UT83v6_CGDV2013.shp"
spatial_reference = 2961

# Approx. coordinates of the Sydney Street Pub
x = 281247.02
y = 4944920.76

radius = 250 # metres

roads = gpd.read_file(test_data)

point = Point(x, y)

buffer = point.buffer(radius)

clipped_roads = gpd.clip(roads, buffer)

# calculate distance for all lines
distances = clipped_roads["geometry"].map(lambda line: point.distance(line))
clipped_roads = clipped_roads.assign(Distance=distances)

# get 4 roads closest to the point
clipped_roads = clipped_roads.sort_values(by="Distance")
closest_roads = clipped_roads.head(4)

Result

print(closest_roads)

Closest Roads


Plotting the result

# plot result
figure, ax = plt.subplots(figsize = (15,7))

# there is probably a more elegant and better way to plot a polygon 
# without converting to a geodataframe first... ;)
# if so, use it as I am not an expert using matplotlib
buffer_gdf = gpd.GeoDataFrame(geometry = [buffer], crs = roads.crs)
buffer_gdf.plot(ax = ax, color = "white", edgecolor = "black")

clipped_roads.plot(ax = ax, color = "black")
closest_roads.plot(ax = ax, color = "red")

point_gdf = gpd.GeoDataFrame(geometry = [point], crs = roads.crs)
point_gdf.plot(ax = ax, color = "blue")

plt.savefig("result.png")

Result


Projecting Coordinates

Following an example how to project coordinates:

from shapely.ops import transform
from shapely.geometry import Point

x = -65.75754
y = 44.624347

wgs84 = pyproj.CRS("EPSG:4326")
utm = pyproj.CRS("EPSG:2961")

project = pyproj.Transformer.from_crs(wgs84, utm, always_xy=True).transform
point = transform(project, Point(x, y))

See also:

1
  • Thank you so much for the help. I think this will work. You are right about the coordinates. but the above code snippet I shared using ArcPy works for me. but I was looking for something with open source. Thanks a lot.
    – alauddin
    Mar 24 at 7:47

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.