I have created 2 plots as pictured:
The blue (with red outline) polygon should be a 10KM circle around the central dark green point in both plots.
This is correct for the right hand side plot (the purple pin on the red line with a light green circle around it is at 10.02KM).
However this is not the case for the left plot. This polygon was created using Geopandas buffer(). Here the same purple pin at 10.02KM is far to the left of its expected location (again surrounded by a light green circle.) The red outline of this polygon is at ~6.2/6.3KM instead of 10KM as expected.
I need to perform a number of geospatial operations on this polygon data including estimating the population of this polygon, hence I need to use GeoPandas.
It appears that my current use of buffer() does not produce a 10,000 Metre buffer radius as expected, but rather one of ~6200-6300 Metres.
How can I make my Geopandas code produce the desired polygon, as illustrated on the right plot?
My code:
# CREATE DF WITH SINGLE ROW
df = pd.DataFrame([[51.502687, -3.538329, 2242, 1, 47, 10.00]],
columns = ["lat","long","num_of_trucks","performance","num_of_routes","avg_agent_dist_trav_km"]
)
# CONVERT ABOVE DF TO GEODF USING LONGITUDE AND LATITUDE:
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["long"],df["lat"]),
crs={"init": "epsg:4326"})
# REPROJECT TO CRS THAT USES METRES SO CAN CALCULATE BUFFER RADIUS
gdf.to_crs(epsg=3395,inplace=True)
# CALCULATE 10KM BUFFER RADIUS <<< THIS STEP SEEMS TO BE CAUSE OF MY ISSUE
gdf["geometry"] = gdf.geometry.buffer(10000)
gdf["current_buffer_radius_in_metres"] = 10000
# CHANGE THE CRS BACK TO ALLOW PLOTTING IN GEOVIEWS
gdf.to_crs(epsg=4326,inplace=True)
# PRODUCE POLYGON FOR LEFT PLOT. THIS IS THE INCORRECT POLYGON:
polyplot = gv.Polygons(gdf).opts(tools=["hover","tap"],line_color="red",line_width=3,fill_alpha=0.35)
# PRODUCE DF TO PRODUCE TRACKER POINT PINS PLOT.
trackersdf = pd.DataFrame([[51.605475, -3.43845, "Portishead Depo Bridgend", 13.370608, 8.30811],
[51.591792, -3.242236, "Portishead Depo Bridgend", 22.806912, 14.171558],
[51.511891, -3.581327, "Portishead Depo Bridgend", 3.156072, 1.961092],
[51.512954, -3.546054, "Portishead Depo Bridgend", 1.261932, 0.784128],
[51.722711, -3.248196, "Portishead Depo Bridgend", 31.672704, 19.680506],
[51.562025, -3.27612, "Portishead Depo Bridgend", 19.355701, 12.027075],
[51.530072, -3.67591, "Portishead Depo Bridgend", 10.024503, 6.228937],
[51.667927, -3.52239, "Portishead Depo Bridgend", 18.417676, 11.444213],
[51.666187, -3.3038, "Portishead Depo Bridgend", 24.395822, 15.158861],
[51.644319, -3.644031, "Portishead Depo Bridgend", 17.378509, 10.798505]],
columns = ["lat", "lon", "depo_name", "tracker_dist_trav_km", "tracker_dist_trav_miles"])
# PRODUCE TRACKER POINT PINS PLOT (USED FOR BOTH PLOTS):
tracker_pins = gv.Points(trackersdf,kdims = ["lon","lat"],
vdims = ["tracker_dist_trav_km","tracker_dist_trav_miles"],
label = "Trackers").opts(tools=["hover","tap"],
color="purple",size = 8,
line_color="white",
hover_color="red",
hover_line_color="white",line_width=0.1,
legend_position="top")
# PRODUCE DEPO POINTS PIN PLOT (USED FOR BOTH PLOTS):
depo_pin = gv.Points(df,
kdims = ["long","lat"],
vdims = ["num_of_trucks","performance","num_of_routes","avg_agent_dist_trav_km"]
).opts(tools=["hover","tap"],
color="green",size = 15,
line_color="white",hover_color="orange",
hover_line_color="white",line_width=2)
### PRODUCE POLYGON FOR THE RIGHT PLOT. THIS PRODUCES THE CORRECT POLYGON, HOWEVER I WISH TO REPRODUCE THIS USING GEOPANDAS.
### THIS IS ADAPTED FROM:
### https://gis.stackexchange.com/questions/289044/creating-buffer-circle-x-kilometers-from-point-using-python
proj_wgs84 = pyproj.Proj('+proj=longlat +datum=WGS84')
def geodesic_point_buffer(index):
lat = df.iloc[index]["lat"]
lon = df.iloc[index]["long"]
km = df.iloc[index]["avg_agent_dist_trav_km"]
# Azimuthal equidistant projection
aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
project = partial(
pyproj.transform,
pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
proj_wgs84)
buf = Point(0, 0).buffer(km * 1000) # DISTANCE IN METRES
# THIS WILL RETURN THE SET OF POINTS (LONG&LAT) THAT TOGETHER MAKE UP THE RADIUS AROUND EACH DEPO
# RETURNED TUPLE 1. DATA:
tracker_radius_coords = transform(project, buf).exterior.coords[:]
# RETURNED TUPLE 2. PLOT:
radius = gv.Points(tracker_radius_coords).opts(color="red",size=5)
# RETURNED TUPLE 3. PLOT:
radius_poly = gv.Polygons(tracker_radius_coords).opts(tools=["hover","tap"],line_color="red",line_width=3,fill_alpha=0.35)
# RETURN MULTIPLE COMPONENTS SO THAT I CAN USE THEM ELSEWHERE
return tracker_radius_coords, radius, radius_poly
# PRODUCE THE RIGHT POLYGON BY CALLING THE ABOVE FUNCTION:
radius_poly = geodesic_point_buffer(0)[2]
# PRODUCE LEFT PLOT. NB THIS IS WRONG.
plot = tracker_pins * depo_pin * polyplot * gvts.OSM
plot.opts(width=800,height=900)
# PRODUCE RIGHT PLOT. NB THIS IS RIGHT, BUT DOES NOT USE GEOPANDAS AS REQUIRED.
plot2 = tracker_pins * depo_pin * radius_poly * gvts.OSM
plot2.opts(width=800,height=900)
# USE A PYVIZ PANEL TO SHOW LEFT AND RIGHT PLOT SIDE BY SIDE AS PER ABOVE PICTURE:
publish = pn.Column(
pn.Row("## Geopandas Version, 10000 metre Buffer Radius used. INCORRECT POLYGON PRODUCED", plot,
"## Non Geopandas Version, 10km Polygon produced. CORRECT POLYGON PRODUCED", plot2))
publish.show(port = 10075,open=True)
