How to Test if Points Intersect a Polygon or Buffer Zone Using GeoPandas and Leaflet

Published: August 24, 2026
Updated: September 02, 2026

Published: August 24, 2026

Updated: September 02, 2026

Introduction

When working with geospatial data in Python, a common task is to test whether points fall inside a polygon.

For example, you may want to know:

  • if a fire detection falls inside a known persistent thermal anomaly area;
  • if a point intersects a solar farm polygon;
  • if a point is outside the polygon but still inside a buffer zone;
  • if random test points intersect a geometry before using the method on a larger dataset.

In this tutorial, we will use GeoPandas, Shapely, and Folium to test whether points intersect a polygon and a buffered polygon.

Table of contents

We will use a polygon extracted from gdf_topaz, representing the largest polygon from a Topaz Solar Farm geometry. Then we will create points around it and add Boolean columns indicating whether each point intersects the polygon or its buffer.

Goal

The goal is to create a GeoDataFrame with points and add columns like this:

1
point_id | longitude | latitude | intersects_polygon | intersects_buffer

Where:

  • intersects_polygon = True means the point intersects the original polygon.
  • intersects_buffer = True means the point intersects the polygon or the surrounding buffer zone.
  • intersects_polygon = False and intersects_buffer = True means the point is outside the polygon but close to it.
  • intersects_buffer = False means the point is outside both the polygon and the buffer.

Required packages

1
2
3
4
5
6
7
8
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import folium

from shapely.geometry import Point, Polygon

If needed, install the packages with:

1
# %pip install geopandas shapely pandas numpy folium mapclassify

Define a polygon from coordinates

The coordinates below are in longitude/latitude order.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
topaz_polygon_coordinates = [
    (-120.028905, 35.360990),
    (-120.027158, 35.360990),
    (-120.027210, 35.353300),
    (-120.032477, 35.353315),
    (-120.032578, 35.359217),
    (-120.037711, 35.359246),
    (-120.037687, 35.368248),
    (-120.039398, 35.368279),
    (-120.039423, 35.378957),
    (-120.037763, 35.378977),
    (-120.037782, 35.380434),
    (-120.027230, 35.380497),
    (-120.027087, 35.363782),
    (-120.028902, 35.363782),
    (-120.028905, 35.360990),
]

Make sure the polygon is closed. A closed polygon means the first and last coordinates are identical.

1
2
if topaz_polygon_coordinates[0] != topaz_polygon_coordinates[-1]:
    topaz_polygon_coordinates.append(topaz_polygon_coordinates[0])

Now create a Shapely polygon:

1
topaz_polygon = Polygon(topaz_polygon_coordinates)

Then create a GeoDataFrame:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
gdf_topaz_polygon = gpd.GeoDataFrame(
    {
        "name": ["Topaz Solar Farm - extracted largest polygon"],
        "source": ["Largest polygon extracted from gdf_topaz"],
    },
    geometry=[topaz_polygon],
    crs="EPSG:4326",
)

gdf_topaz_polygon
name source geometry area_km2
Topaz Solar Farm - extracted largest polygon Largest polygon extracted from gdf_topaz POLYGON ((-120.028905 35.36099, -120.027158 35.36099, -120.02721 35.3533, -120.032477 35.353315, -120.032578 35.359217, -120.037711 35.359246, -120.037687 35.368248, -120.039398 35.368279, -120.039423 35.378957, -120.037763 35.378977, -120.037782 35.380434, -120.02723 35.380497, -120.027087 35.363782, -120.028902 35.363782, -120.028905 35.36099)) 2.7147696427931787

The CRS is EPSG:4326, meaning the coordinates are longitude and latitude.

Check the polygon coordinates

You can extract the polygon exterior coordinates back into a pandas DataFrame:

1
2
3
4
5
6
coords_df = pd.DataFrame(
    list(gdf_topaz_polygon.geometry.iloc[0].exterior.coords),
    columns=["longitude", "latitude"]
)

coords_df
longitude latitude
-120.028905 35.36099
-120.027158 35.36099
-120.02721 35.3533
-120.032477 35.353315
-120.032578 35.359217
-120.037711 35.359246
-120.037687 35.368248
-120.039398 35.368279
-120.039423 35.378957
-120.037763 35.378977
-120.037782 35.380434
-120.02723 35.380497
-120.027087 35.363782
-120.028902 35.363782
-120.028905 35.36099

This is useful if you want to inspect or export the polygon coordinates.

Compute the polygon area in square kilometers

Do not calculate area directly in latitude/longitude coordinates. In EPSG:4326, the geometry units are degrees, not meters.

Instead, estimate a local projected CRS and calculate the area there.

1
2
3
utm_crs = gdf_topaz_polygon.estimate_utm_crs()

print("Estimated UTM CRS:", utm_crs)

Project the polygon to the estimated UTM CRS:

1
gdf_topaz_polygon_m = gdf_topaz_polygon.to_crs(utm_crs)

Now calculate the area in square kilometers:

1
2
3
4
5
gdf_topaz_polygon["area_km2"] = (
    gdf_topaz_polygon_m.geometry.area.values / 1_000_000
)

gdf_topaz_polygon[["name", "area_km2", "geometry"]]

Estimated UTM CRS: EPSG:32610

name area_km2 geometry
Topaz Solar Farm - extracted largest polygon 2.7147696427931787 POLYGON ((-120.028905 35.36099, -120.027158 35.36099, -120.02721 35.3533, -120.032477 35.353315, -120.032578 35.359217, -120.037711 35.359246, -120.037687 35.368248, -120.039398 35.368279, -120.039423 35.378957, -120.037763 35.378977, -120.037782 35.380434, -120.02723 35.380497, -120.027087 35.363782, -120.028902 35.363782, -120.028905 35.36099))

Create manual test points

Before generating many random points, it is useful to create a few manual points.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
df_test_points = pd.DataFrame({
    "point_id": [1, 2, 3, 4, 5, 6],
    "description": [
        "likely inside polygon",
        "inside or close to polygon",
        "east of polygon",
        "west of polygon",
        "north of polygon",
        "south of polygon",
    ],
    "longitude": [
        -120.0320,
        -120.0285,
        -120.0245,
        -120.0420,
        -120.0330,
        -120.0330,
    ],
    "latitude": [
        35.3650,
        35.3700,
        35.3660,
        35.3680,
        35.3830,
        35.3510,
    ],
})

Convert the pandas DataFrame to a GeoDataFrame:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
gdf_test_points = gpd.GeoDataFrame(
    df_test_points,
    geometry=gpd.points_from_xy(
        df_test_points["longitude"],
        df_test_points["latitude"],
    ),
    crs="EPSG:4326",
)

gdf_test_points
point_id description longitude latitude geometry intersects_polygon intersects_buffer intersection_class
1 likely inside polygon -120.032 35.365 POINT (-120.032 35.365) True True Inside polygon
2 inside or close to polygon -120.0285 35.37 POINT (-120.0285 35.37) True True Inside polygon
3 east of polygon -120.0245 35.366 POINT (-120.0245 35.366) False True Inside buffer only
4 west of polygon -120.042 35.368 POINT (-120.042 35.368) False True Inside buffer only
5 north of polygon -120.033 35.383 POINT (-120.033 35.383) False True Inside buffer only
6 south of polygon -120.033 35.351 POINT (-120.033 35.351) False True Inside buffer only

Test whether points intersect the polygon

Extract the polygon geometry:

1
topaz_geom = gdf_topaz_polygon.geometry.iloc[0]

Then use intersects():

1
2
3
gdf_test_points["intersects_polygon"] = (
    gdf_test_points.geometry.intersects(topaz_geom)
)

Display the result:

1
2
3
4
5
6
7
8
9
gdf_test_points[
    [
        "point_id",
        "description",
        "longitude",
        "latitude",
        "intersects_polygon",
    ]
]
point_id description longitude latitude intersects_polygon
1 likely inside polygon -120.032 35.365 True
2 inside or close to polygon -120.0285 35.37 True
3 east of polygon -120.0245 35.366 False
4 west of polygon -120.042 35.368 False
5 north of polygon -120.033 35.383 False
6 south of polygon -120.033 35.351 False

The new column intersects_polygon contains True or False.

Why use intersects?

For points and polygons, you can use several spatial predicates:

1
2
3
gdf_test_points.geometry.within(topaz_geom)
gdf_test_points.geometry.intersects(topaz_geom)
gdf_test_points.geometry.touches(topaz_geom)

The difference matters:

  • within is true only when the point is strictly inside the polygon.
  • touches is true when the point is exactly on the polygon boundary.
  • intersects is true if the point is inside or touches the polygon.

For many Earth observation applications, intersects is a good choice because it includes boundary cases.

Create a buffer around the polygon

A buffer is useful when the polygon boundary is approximate, or when the point location has some uncertainty.

For example, a satellite fire pixel, a persistent thermal anomaly, or a geolocation point may be slightly shifted. A buffer allows us to include points close to the polygon.

Important: create the buffer in a projected CRS using meters.

1
buffer_m = 1000

Project the polygon and points:

1
2
gdf_topaz_polygon_m = gdf_topaz_polygon.to_crs(utm_crs)
gdf_test_points_m = gdf_test_points.to_crs(utm_crs)

Create the buffer:

1
topaz_buffer_m = gdf_topaz_polygon_m.geometry.iloc[0].buffer(buffer_m)

Create a GeoDataFrame for the buffer:

1
2
3
4
5
gdf_topaz_buffer_m = gpd.GeoDataFrame(
    {"name": [f"Topaz Solar Farm {buffer_m} m buffer"]},
    geometry=[topaz_buffer_m],
    crs=utm_crs,
)

Convert the buffer back to latitude/longitude for mapping:

1
gdf_topaz_buffer = gdf_topaz_buffer_m.to_crs("EPSG:4326")

Now test whether each point intersects the buffer:

1
2
3
4
5
6
7
gdf_test_points_m["intersects_buffer"] = (
    gdf_test_points_m.geometry.intersects(topaz_buffer_m)
)

gdf_test_points["intersects_buffer"] = (
    gdf_test_points_m["intersects_buffer"].values
)

Display the result:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
gdf_test_points[
    [
        "point_id",
        "description",
        "longitude",
        "latitude",
        "intersects_polygon",
        "intersects_buffer",
    ]
]
point_id description longitude latitude intersects_polygon intersects_buffer
1 likely inside polygon -120.032 35.365 True True
2 inside or close to polygon -120.0285 35.37 True True
3 east of polygon -120.0245 35.366 False True
4 west of polygon -120.042 35.368 False True
5 north of polygon -120.033 35.383 False True
6 south of polygon -120.033 35.351 False True

Add a readable class column

The Boolean columns are useful, but a readable class is better for mapping.

1
2
3
4
5
6
7
def classify_intersection(row):
    if row["intersects_polygon"]:
        return "Inside polygon"
    elif row["intersects_buffer"]:
        return "Inside buffer only"
    else:
        return "Outside"

Apply the function:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
gdf_test_points["intersection_class"] = gdf_test_points.apply(
    classify_intersection,
    axis=1,
)

gdf_test_points[
    [
        "point_id",
        "description",
        "intersects_polygon",
        "intersects_buffer",
        "intersection_class",
    ]
]
point_id description intersects_polygon intersects_buffer intersection_class
1 likely inside polygon True True Inside polygon
2 inside or close to polygon True True Inside polygon
3 east of polygon False True Inside buffer only
4 west of polygon False True Inside buffer only
5 north of polygon False True Inside buffer only
6 south of polygon False True Inside buffer only

Generate random points around the polygon

Now let’s generate random points around the polygon.

First, get the polygon bounds:

1
min_lon, min_lat, max_lon, max_lat = gdf_topaz_polygon.total_bounds

Add some padding around the polygon:

1
2
padding_lon = 0.015
padding_lat = 0.015

Generate random points:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
rng = np.random.default_rng(42)

n_points = 100

random_lons = rng.uniform(
    min_lon - padding_lon,
    max_lon + padding_lon,
    n_points,
)

random_lats = rng.uniform(
    min_lat - padding_lat,
    max_lat + padding_lat,
    n_points,
)

Create a DataFrame:

1
2
3
4
5
df_random_points = pd.DataFrame({
    "point_id": range(1, n_points + 1),
    "longitude": random_lons,
    "latitude": random_lats,
})

Convert to a GeoDataFrame:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
gdf_random_points = gpd.GeoDataFrame(
    df_random_points,
    geometry=gpd.points_from_xy(
        df_random_points["longitude"],
        df_random_points["latitude"],
    ),
    crs="EPSG:4326",
)

gdf_random_points.head()
point_id longitude latitude geometry intersects_polygon intersects_buffer intersection_class
1 -120.02165679672834 35.3902680897664 POINT (-120.02165679672834 35.3902680897664) False False Outside
2 -120.03584264237466 35.37832114893257 POINT (-120.03584264237466 35.37832114893257) True True Inside polygon
3 -120.01807339846263 35.3535069641856 POINT (-120.01807339846263 35.3535069641856) False True Inside buffer only
4 -120.02489922712174 35.39373398125516 POINT (-120.02489922712174 35.39373398125516) False False Outside
5 -120.05043590779982 35.38284221545413 POINT (-120.05043590779982 35.38284221545413) False False Outside

Test random points against the polygon and buffer

Test the original polygon:

1
2
3
gdf_random_points["intersects_polygon"] = (
    gdf_random_points.geometry.intersects(topaz_geom)
)

Test the buffer:

1
2
3
4
5
6
7
8
9
gdf_random_points_m = gdf_random_points.to_crs(utm_crs)

gdf_random_points_m["intersects_buffer"] = (
    gdf_random_points_m.geometry.intersects(topaz_buffer_m)
)

gdf_random_points["intersects_buffer"] = (
    gdf_random_points_m["intersects_buffer"].values
)

Add the readable class:

1
2
3
4
5
6
gdf_random_points["intersection_class"] = gdf_random_points.apply(
    classify_intersection,
    axis=1,
)

gdf_random_points.head(15)
point_id longitude latitude geometry intersects_polygon intersects_buffer intersection_class
1 -120.02165679672834 35.3902680897664 POINT (-120.02165679672834 35.3902680897664) False False Outside
2 -120.03584264237466 35.37832114893257 POINT (-120.03584264237466 35.37832114893257) True True Inside polygon
3 -120.01807339846263 35.3535069641856 POINT (-120.01807339846263 35.3535069641856) False True Inside buffer only
4 -120.02489922712174 35.39373398125516 POINT (-120.02489922712174 35.39373398125516) False False Outside
5 -120.05043590779982 35.38284221545413 POINT (-120.05043590779982 35.38284221545413) False False Outside
6 -120.0131190521211 35.37930396814932 POINT (-120.0131190521211 35.37930396814932) False False Outside
7 -120.02219938957653 35.36400212983811 POINT (-120.02219938957653 35.36400212983811) False True Inside buffer only
8 -120.02114418157178 35.35387140061285 POINT (-120.02114418157178 35.35387140061285) False True Inside buffer only
9 -120.04899918124704 35.34381327386229 POINT (-120.04899918124704 35.34381327386229) False False Outside
10 -120.03535546093325 35.38992614927512 POINT (-120.03535546093325 35.38992614927512) False False Outside
11 -120.03872489484608 35.36436903644961 POINT (-120.03872489484608 35.36436903644961) False True Inside buffer only
12 -120.0151874774321 35.34987457737619 POINT (-120.0151874774321 35.34987457737619) False False Outside
13 -120.02716432627626 35.355799801031544 POINT (-120.02716432627626 35.355799801031544) False True Inside buffer only
14 -120.01959056434056 35.37142962168477 POINT (-120.01959056434056 35.37142962168477) False True Inside buffer only
15 -120.03565061647845 35.34841087286577 POINT (-120.03565061647845 35.34841087286577) False True Inside buffer only

Count points by class

1
2
3
4
5
6
7
8
9
summary = (
    gdf_random_points
    .groupby("intersection_class")
    .size()
    .reset_index(name="n_points")
    .sort_values("intersection_class")
)

summary
intersection_class n_points
Inside buffer only 53
Inside polygon 11
Outside 36

This gives a quick summary of how many points are:

  • inside the polygon;
  • inside the buffer only;
  • outside both.

Spatial join alternative

For a single polygon, this is enough:

1
gdf_random_points.geometry.intersects(topaz_geom)

But if you have many polygons, gpd.sjoin() is often more convenient.

1
2
3
4
5
6
matches = gpd.sjoin(
    gdf_random_points,
    gdf_topaz_polygon[["name", "geometry"]],
    how="left",
    predicate="intersects",
)

Inspect the result:

1
2
3
4
5
6
7
8
9
matches[
    [
        "point_id",
        "longitude",
        "latitude",
        "intersects_polygon",
        "name",
    ]
].head(20)
point_id longitude latitude intersects_polygon name
1 -120.02165679672834 35.3902680897664 False nan
2 -120.03584264237466 35.37832114893257 True Topaz Solar Farm - extracted largest polygon
3 -120.01807339846263 35.3535069641856 False nan
4 -120.02489922712174 35.39373398125516 False nan
5 -120.05043590779982 35.38284221545413 False nan
6 -120.0131190521211 35.37930396814932 False nan
7 -120.02219938957653 35.36400212983811 False nan
8 -120.02114418157178 35.35387140061285 False nan
9 -120.04899918124704 35.34381327386229 False nan
10 -120.03535546093325 35.38992614927512 False nan
11 -120.03872489484608 35.36436903644961 False nan
12 -120.0151874774321 35.34987457737619 False nan
13 -120.02716432627626 35.355799801031544 False nan
14 -120.01959056434056 35.37142962168477 False nan
15 -120.03565061647845 35.34841087286577 False nan
16 -120.04480262147452 35.38729576720723 False nan
17 -120.0309440984569 35.38168504154798 False nan
18 -120.05172123264558 35.37945112269152 False nan
19 -120.01938440670251 35.36301442559601 False nan
20 -120.02768085599877 35.37418018376166 True Topaz Solar Farm - extracted largest polygon

Points intersecting the polygon will receive the polygon name.

Points outside the polygon will have missing values in the joined polygon columns.

Export polygon, buffer, and points to GeoJSON

GeoJSON is useful because it can be displayed directly in Leaflet.

1
2
output_dir = Path("topaz_polygon_buffer_leaflet_example")
output_dir.mkdir(exist_ok=True)

Save the polygon:

1
2
3
4
gdf_topaz_polygon.to_file(
    output_dir / "topaz_solar_farm_polygon.geojson",
    driver="GeoJSON",
)

Save the buffer:

1
2
3
4
gdf_topaz_buffer.to_file(
    output_dir / "topaz_solar_farm_buffer.geojson",
    driver="GeoJSON",
)

Save the random points:

1
2
3
4
gdf_random_points.to_file(
    output_dir / "random_points_topaz.geojson",
    driver="GeoJSON",
)

Save the manual test points:

1
2
3
4
gdf_test_points.to_file(
    output_dir / "manual_test_points_topaz.geojson",
    driver="GeoJSON",
)

Check the files:

1
2
3
print("Created files:")
for path in sorted(output_dir.glob("*.geojson")):
    print(path)

Create a Leaflet map with Folium

Folium allows us to create a Leaflet map directly from Python.

centroid = gdf_topaz_polygon.geometry.iloc[0].centroid

1
2
3
4
5
m = folium.Map(
    location=[centroid.y, centroid.x],
    zoom_start=14,
    tiles="OpenStreetMap",
)

Add the buffer first:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
folium.GeoJson(
    gdf_topaz_buffer,
    name=f"{buffer_m} m buffer",
    style_function=lambda feature: {
        "color": "#9333ea",
        "weight": 2,
        "fillColor": "#9333ea",
        "fillOpacity": 0.15,
    },
).add_to(m)

Add the original polygon:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
folium.GeoJson(
    gdf_topaz_polygon,
    name="Topaz polygon",
    style_function=lambda feature: {
        "color": "#2563eb",
        "weight": 3,
        "fillColor": "#2563eb",
        "fillOpacity": 0.25,
    },
    tooltip=folium.GeoJsonTooltip(fields=["name", "area_km2"]),
).add_to(m)

Define colors for the points:

1
2
3
4
5
color_lookup = {
    "Inside polygon": "green",
    "Inside buffer only": "orange",
    "Outside": "red",
}

Add the random points:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
for _, row in gdf_random_points.iterrows():
    color = color_lookup[row["intersection_class"]]

    popup_html = (
        f"<b>Random point ID:</b> {row['point_id']}<br>"
        f"<b>Intersects polygon:</b> {row['intersects_polygon']}<br>"
        f"<b>Intersects buffer:</b> {row['intersects_buffer']}<br>"
        f"<b>Class:</b> {row['intersection_class']}"
    )

    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],
        radius=4,
        color=color,
        fill=True,
        fill_color=color,
        fill_opacity=0.85,
        popup=folium.Popup(popup_html, max_width=260),
    ).add_to(m)

Add manual test points with larger markers:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
for _, row in gdf_test_points.iterrows():
    color = color_lookup[row["intersection_class"]]

    popup_html = (
        f"<b>Manual point ID:</b> {row['point_id']}<br>"
        f"<b>Description:</b> {row['description']}<br>"
        f"<b>Intersects polygon:</b> {row['intersects_polygon']}<br>"
        f"<b>Intersects buffer:</b> {row['intersects_buffer']}<br>"
        f"<b>Class:</b> {row['intersection_class']}"
    )

    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],
        radius=8,
        color="black",
        weight=2,
        fill=True,
        fill_color=color,
        fill_opacity=0.95,
        popup=folium.Popup(popup_html, max_width=280),
    ).add_to(m)

Add a layer control:

1
2
3
folium.LayerControl().add_to(m)

m

How to Test if Points Intersect a Polygon or Buffer Zone Using GeoPandas and Leaflet
How to Test if Points Intersect a Polygon or Buffer Zone Using GeoPandas and Leaflet

The map colors are:

  • green: point intersects the polygon;
  • orange: point intersects only the buffer;
  • red: point is outside both.

Save the map as HTML

1
2
3
4
5
html_map_path = output_dir / "topaz_polygon_buffer_intersection_map.html"

m.save(html_map_path)

print("Saved map to:", html_map_path)

You can open the saved HTML file in a browser.

Create a standalone Leaflet page

If you want to create a pure Leaflet page, you can use the exported GeoJSON files.

Create a file called:

1
index.html

inside the folder:

1
topaz_polygon_buffer_leaflet_example

Then use this HTML template:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="UTF-8" />
  <title>Topaz Solar Farm Polygon Buffer Example</title>
  <meta name="viewport" content="width=device-width, initial-scale=1.0" />

  <link
    rel="stylesheet"
    href="https://unpkg.com/leaflet@1.9.4/dist/leaflet.css"
  />

  <style>
    body {
      margin: 0;
      font-family: Arial, sans-serif;
      background: #f4f4f4;
    }

    header {
      padding: 16px 24px;
      background: #111827;
      color: white;
    }

    header h1 {
      margin: 0;
      font-size: 22px;
    }

    header p {
      margin: 6px 0 0 0;
      color: #d1d5db;
      font-size: 14px;
    }

    #map {
      width: 100%;
      height: 720px;
    }

    .legend {
      position: absolute;
      right: 16px;
      bottom: 24px;
      z-index: 1000;
      background: white;
      padding: 12px 14px;
      border-radius: 8px;
      box-shadow: 0 2px 8px rgba(0,0,0,0.25);
      font-size: 13px;
      line-height: 1.5;
    }

    .legend-item {
      display: flex;
      align-items: center;
      gap: 8px;
    }

    .dot {
      width: 12px;
      height: 12px;
      border-radius: 50%;
      display: inline-block;
    }

    .dot-inside {
      background: #16a34a;
    }

    .dot-buffer {
      background: #f59e0b;
    }

    .dot-outside {
      background: #dc2626;
    }

    .line-polygon {
      width: 18px;
      height: 4px;
      background: #2563eb;
      display: inline-block;
    }

    .line-buffer {
      width: 18px;
      height: 4px;
      background: #9333ea;
      display: inline-block;
      opacity: 0.5;
    }
  </style>
</head>

<body>

<header>
  <h1>Testing Point Intersections with a Topaz Solar Farm Polygon</h1>
  <p>
    Green points intersect the polygon. Orange points intersect only the buffer.
    Red points are outside both.
  </p>
</header>

<div id="map"></div>

<div class="legend">
  <div class="legend-item">
    <span class="dot dot-inside"></span>
    Point intersects polygon
  </div>
  <div class="legend-item">
    <span class="dot dot-buffer"></span>
    Point intersects buffer only
  </div>
  <div class="legend-item">
    <span class="dot dot-outside"></span>
    Point outside polygon and buffer
  </div>
  <div class="legend-item">
    <span class="line-polygon"></span>
    Topaz polygon
  </div>
  <div class="legend-item">
    <span class="line-buffer"></span>
    Buffer zone
  </div>
</div>

<script src="https://unpkg.com/leaflet@1.9.4/dist/leaflet.js"></script>

<script>
  const map = L.map("map").setView([35.3669, -120.0333], 14);

  L.tileLayer("https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png", {
    maxZoom: 19,
    attribution: "&copy; OpenStreetMap contributors"
  }).addTo(map);

  function pointColor(props) {
    if (props.intersects_polygon === true || props.intersects_polygon === "True") {
      return "#16a34a";
    }

    if (props.intersects_buffer === true || props.intersects_buffer === "True") {
      return "#f59e0b";
    }

    return "#dc2626";
  }

  fetch("topaz_solar_farm_buffer.geojson")
    .then(response => response.json())
    .then(data => {
      L.geoJSON(data, {
        style: {
          color: "#9333ea",
          weight: 2,
          fillColor: "#9333ea",
          fillOpacity: 0.15
        }
      }).addTo(map);
    });

  fetch("topaz_solar_farm_polygon.geojson")
    .then(response => response.json())
    .then(data => {
      const polygonLayer = L.geoJSON(data, {
        style: {
          color: "#2563eb",
          weight: 3,
          fillColor: "#2563eb",
          fillOpacity: 0.25
        },
        onEachFeature: function(feature, layer) {
          layer.bindPopup(feature.properties.name);
        }
      }).addTo(map);

      map.fitBounds(polygonLayer.getBounds(), {
        padding: [30, 30]
      });
    });

  fetch("random_points_topaz.geojson")
    .then(response => response.json())
    .then(data => {
      L.geoJSON(data, {
        pointToLayer: function(feature, latlng) {
          const color = pointColor(feature.properties);

          return L.circleMarker(latlng, {
            radius: 5,
            color: color,
            fillColor: color,
            fillOpacity: 0.85,
            weight: 1
          });
        },
        onEachFeature: function(feature, layer) {
          const props = feature.properties;

          const popupContent = `
            <strong>Point ID:</strong> ${props.point_id}<br>
            <strong>Intersects polygon:</strong> ${props.intersects_polygon}<br>
            <strong>Intersects buffer:</strong> ${props.intersects_buffer}<br>
            <strong>Class:</strong> ${props.intersection_class}
          `;

          layer.bindPopup(popupContent);
        }
      }).addTo(map);
    });
</script>

</body>
</html>

Run the Leaflet example locally

Because the HTML file uses fetch() to load local GeoJSON files, it is better to use a local web server.

From the terminal:

1
2
cd topaz_polygon_buffer_leaflet_example
python -m http.server 8000

Then open:

1
http://localhost:8000/index.html

Common mistake: creating buffers in latitude/longitude

This is a common mistake:

1
gdf_topaz_polygon.geometry.buffer(1000)

If the CRS is EPSG:4326, the buffer distance is interpreted in degrees, not meters.

The safer workflow is:

1
2
3
4
5
utm_crs = gdf_topaz_polygon.estimate_utm_crs()

gdf_topaz_polygon_m = gdf_topaz_polygon.to_crs(utm_crs)

buffer_polygon_m = gdf_topaz_polygon_m.geometry.iloc[0].buffer(1000)

Then convert the buffer back to longitude/latitude for Leaflet:

1
2
3
4
5
gdf_topaz_buffer = gpd.GeoDataFrame(
    {"name": ["1000 m buffer"]},
    geometry=[buffer_polygon_m],
    crs=utm_crs,
).to_crs("EPSG:4326")

Conclusion

GeoPandas makes it easy to test whether points intersect a polygon or a buffer zone.

The main workflow is:

  1. Create or load a polygon.
  2. Convert point coordinates to a GeoDataFrame.
  3. Use geometry.intersects(polygon) to create a Boolean column.
  4. Project to a meter-based CRS before creating buffers.
  5. Use the buffered polygon to create a second Boolean column.
  6. Export the result to GeoJSON.
  7. Display the result with Folium or Leaflet.

This approach is useful for many Earth observation workflows, including persistent thermal anomaly analysis, fire detection validation, solar farm mapping, and spatial filtering of satellite observations.

References

Links Site
GeoPandas documentation GeoPandas
GeoPandas spatial joins GeoPandas
GeoPandas estimate_utm_crs GeoPandas
Shapely Polygon documentation Shapely
Shapely intersects documentation Shapely
Folium documentation Folium
Leaflet GeoJSON tutorial Leaflet
OpenStreetMap OpenStreetMap
Author profile-image

Benjamin

Earth observation scientist and developer, I explore satellite data, fire detection, and environmental analytics using Python and machine learning. Passionate about open science and sharing practical tools.

Interests
  • Open Science
  • Remote sensing
  • Python
  • Machine Learning
  • Teaching
Jupyter Notebook

Run and modify the Python examples from this article.

Available to Moonbooks members
Creative Commons License

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Image

of