2

My question is similar to this one, but I am having trouble applying the solution to my case.

I have two data sets: one contains a bunch of latitude / longitude points that I've converted into a Point geometry column. The other is a shapefile containing information on primary roads in the US from the Census website. My goal is to compute the distance between each point and the closest primary road within some buffer area.

I think there are two sources of confusion for me:

  1. I'm not sure if I'm using the correct CRS or if I need to be projecting my latitude/longitude points. I find projection and when to use it a bit confusing still.
  2. Am I able to create a buffer and clip the roads around each point separately?

Right now what is happening is that when I try to clip my computer terminates the process after a while. Trying on individual points I get null results.

Here is what I'm trying right now:

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
from shapely import geometry
from shapely.ops import transform
from shapely.geometry import Point


# read in my points data
df = pd.read_csv(myfile)
geo_df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
geo_df = geo_df.set_crs('EPSG:3857')

# read in the roads data
primary_roads = gpd.read_file(tl_2021_us_primaryroads.shp)
primary_roads = primary_roads.to_crs('EPSG:3857')

# define a buffer radius
radius = 1000

# get buffer around each point
geo_df['buffer'] = [point.buffer(radius) for point in geo_df['geometry']]

# create new geometry of clipped roads for each point
geo_df['clipped_roads'] = [gpd.clip(primary_roads, buffer) for buffer in geo_df['buffer']]

The last line is where I'm failing to make progress.

Here is a sample data if it helps:

    latitude  longitude                geometry
0  33.516935 -86.779229  POINT (-86.779 33.517)
1  32.415983 -86.295076  POINT (-86.295 32.416)
2  34.414842 -87.045289  POINT (-87.045 34.415)
3  32.300587 -86.208014  POINT (-86.208 32.301)
4  30.673394 -88.140719  POINT (-88.141 30.673)

The primary roads file can be downloaded here

Maybe I am misunderstanding how clip works.

UPDATE: After some tinkering around, I got the accepted answer to work for me. I ended up using 'ESRI:102003' as my CRS, and this seems to work fine for my purposes. One thing that tripped me up for a bit was that I was setting my points crs to 102003 immediately after reading in the data, which was giving me the wrong coordinates. I needed to first set crs to 'EPSG:4269' and then convert to the other crs.

2
  • 2
    EPSG:3857 is intended to be used for web mapping and it is not good for distances. It is better to use a projected CRS such as UTM which is suitable for such analysis. You need to know the correct zone for your data before projecting the data which depends on the location of your study area.
    – ahmadhanb
    Dec 2, 2022 at 5:52
  • Is there one UTM code that I could use for the entire US area?
    – Sergei
    Dec 2, 2022 at 13:17

1 Answer 1

1

Reproject to a coordinate system of your choice. The example coordinates you are showing are in Alabama US so I choose 32613, I dont know which is best if you have data all over the US.

Buffer, intersect buffers with roads, for each point id find the closest road, merge(=join) the result back to the points:

import geopandas as gpd
import numpy as np

roads = gpd.read_file(r"/home/bera/Downloads/tl_2022_us_primaryroads/tl_2022_us_primaryroads.shp")
roads = roads.to_crs(32616)
roads = roads.drop_duplicates(subset="geometry") #There are some duplicate roads on top of eachother, drop them.

points = gpd.read_file(r"/home/bera/Downloads/tl_2022_us_primaryroads/random_points.shp")
points = points.to_crs(32616)
points["pointid"] = np.arange(points.shape[0]) #Create an id column
points["pointgeometry"] = points["geometry"] #Save the point geometry
points["geometry"] = points.buffer(10000)

#Intersect to clip all roads with each buffer and get all road attributes
clipped = points.overlay(roads, how="intersection", keep_geom_type=False)

closest = []
for group, frame in clipped.groupby("pointid"): #For each pointid and all lines intersecting its buffer
    #Calculate the distances to all roads intersecting the point buffer
    frame["distance"] = frame["pointgeometry"].distance(frame["geometry"]) #That is why we saved the point geometry
    m = frame[frame["distance"]==frame["distance"].min()] #Find the road closest road
    closest.append(m) #m is a dataframe, append it to closest list

closest = gpd.pd.concat(closest) #Create a dataframe from the list

#There are a few points that are exactly the same distance two different roads, I drop them. I dont know what you want to do
closest = closest.drop_duplicates(subset="pointid")

#You can then join the result back to the points
points = points.merge(right=closest[["pointid","LINEARID","distance"]], how="left", on="pointid")

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.