3

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?

enter image description here

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?

3
  • 3
    Geopandas 0.7.0 has clip module designed to do exactly this: geopandas.readthedocs.io/en/latest/gallery/plot_clip.html Commented Feb 24, 2020 at 17:05
  • Thanks! However it gives me Empty GeoDataFrame when I'm clipping Linestrings with Polygon
    – Vilq
    Commented Feb 25, 2020 at 15:39
  • It should not. Do you have the same CRS in both GeoDataFrames? Commented Feb 25, 2020 at 17:25

1 Answer 1

3

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:

enter image description here

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