3

I would like to create an accurate buffer of 5 miles around a coordinate, my current code is: 

cpr_gdf['p_buffer']=cpr_gdf['coordinates'].buffer(5*(1/60))

The coordinates column was created with this code: 

cpr_df['coordinates']=list(zip(cpr_df.sample_longitude_decimal,cpr_df.sample_latitude_decimal))
cpr_df['coordinates']=cpr_df['coordinates'].apply(Point)
cpr_gdf=gpd.GeoDataFrame(cpr_df,geometry='coordinates',crs={'init' :'epsg:4326'})

Thanks for any help!

iacob
  • 20,084
  • 6
  • 92
  • 119
lourew
  • 65
  • 1
  • 1
  • 7
  • Can you clarify what exactly is your question? What about the code example you show does not work (as intended)? – joris Jul 20 '18 at 22:49
  • One remark based on the code: your data is in longitude/latitude. If you want an accurate buffer in miles (meters), you will probably need to convert your data to another coordinate reference system that works in meters. – joris Jul 20 '18 at 22:51
  • Hi, the code has been adjusted so it works :) should I post it in here or just delete the question. – lourew Jul 22 '18 at 09:38
  • You can answer your own question. – joris Jul 22 '18 at 11:11

3 Answers3

4

You need to convert to an equal area projection that is most accurate to where your buffer will be (good resource at https://epsg.io/)

For example, I'm making maps in Michigan, so I'm using EPSG:3174 (which I believe is in meters, correct me if wrong). Given you've already converted your dataframe to a GeoPandas dataframe, you can convert your current projection to 3174 and then create your buffer (converting miles to meters)

cpr_gdf= cpr_gdf.to_crs({'init': 'epsg:3174'})  
buffer_length_in_meters = (5 * 1000) * 1.60934
cpr_gdf['geometry'] = cpr_gdf.geometry.buffer(buffer_length_in_meters)
Geoff Perrin
  • 444
  • 1
  • 6
  • 14
0

You can calculate buffer over points without converting to any other CRS using the function bellow. But it calculates in meters, so if you want to use miles just multiply distance on 1609.34 Here is an example

    from geographiclib.geodesic import Geodesic
    import numpy as np
    from shapely.geometry import Polygon
    import pandas as pd
    import geopandas as gpd
    
    
    def geod_buffer(gdf, distance, resolution=16, geod = Geodesic.WGS84):
        """
        gdf - GeoDataFrame with geometry column
        distance - The radius of the buffer in meters
        resolution – The resolution of the buffer around each vertex
        geod - Define an ellipsoid
        """
        buffer = list()
        for index, row in gdf.iterrows():
            lon1, lat1 = row['geometry'].x, row['geometry'].y
            buffer_ = list()
            for azi1 in np.arange(0, 360, 90/resolution):
                properties = geod.Direct(lat1, lon1, azi1, distance)
                buffer_.append([properties['lon2'], properties['lat2']])
            buffer.append(Polygon(buffer_))
        return buffer
    
    
    locations = pd.DataFrame([
        {
         'longitude': 54.604972,
         'latitude': 18.346815},
        {
         'longitude': 54.605917,
         'latitude': 18.347249}
    ])
    
    locations_gpd = gpd.GeoDataFrame(locations,
                                     geometry=gpd.points_from_xy(locations.longitude, locations.latitude),
                                     crs='epsg:4326').drop(columns=['longitude', 'latitude'])
    
    
    locations_gpd['geometry'] = geod_buffer(locations_gpd, 1000)
Grigory Nevsky
  • 1
  • 1
-1

At the equator, one minute of latitude or longitude is ~ 1.84 km or 1.15 mi (ref).

So if you define your point as P = [y, x] then you can create a buffer around it of lets say 4 minutes which are approximately 5 miles: buffer = 0.04. The bounding box then is easily obtained with

minlat = P[0]-(P[0]*buffer)
maxlat = P[0]+(P[0]*buffer)
minlon = P[1]-(P[1]*buffer)
maxlon = P[1]+(P[1]*buffer)
lorenzori
  • 737
  • 7
  • 23
  • 1
    this only works for certain latitudes, and is an approximation, at best. Results will be wildly skewed at the poles. – Geoff Perrin Jul 26 '19 at 16:46