Python clipping large shapefile by radius [closed]

Question

I cannot find any real answer so I'd like to ask how to crop shapefile with roads in Python.

The goal is to leave roads (nodes) in black circle only, then output that to dataframe in geopandas (possibly save later as shapefile).

Assuming, that I have some place's coordinates in lat + lon and a very big file of roads, how to clip this file by radius, let's say, of 1km?

Theoretically I can create shapefile of the circle basing on coordinates and then crop roads by this shapefile using some external library like ogr2ogr, but it seems to be very slow and I don't want to use subprocess.

Is it possible to do that directly in geopandas? (I can load all roads to pd dataframe first)

Is there any already implemented option to crop roads but leave all these nodes that are even partially in black circle?

Answer 1

Solution is fairly simple, maybe it'll be useful for somebody.

Thanks to @martinfleis and this answer.

First - we have point:

point = Point(-1.7063,55.4127)

Next create a circle basing on point and with buffer - in here 1000 (meters!):

    local_azimuthal_projection = f"+proj=aeqd +R=6371000 +units=m +lat_0={point.y} +lon_0={point.x}"

    wgs84_to_aeqd = partial(
        pyproj.transform,
        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),
        pyproj.Proj(local_azimuthal_projection),
    )

    aeqd_to_wgs84 = partial(
        pyproj.transform,
        pyproj.Proj(local_azimuthal_projection),
        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),
    )

    aeqd_to_27700 = partial(
        pyproj.transform,
        pyproj.Proj(local_azimuthal_projection),
        pyproj.Proj('+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs'),
    )

    point_transformed = transform(wgs84_to_aeqd, point)
    loc_buffer = point_transformed.buffer(1000)

Then apply appropriate transformation:

buffer_wgs84 = transform(aeqd_to_27700, loc_buffer)

load shapefile:

    roadfile = 'files/OSOpenRoads_NU.gml'  # ---> OS is in 27700
    roads   = geopandas.read_file(roadfile)
    clipped = geopandas.clip(roads, buffer_wgs84)

After clipping you can print clipped dataframe or output that.

Tested on Python 3.8 and GeoPandas 0.7 so to use that you'll have to have newer versions.

HOWEVER - sometimes behaviour is a little bit weird, some streets are in the list even if they are outside of boundary, some are not included: