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)