diff --git a/content/algorithms/geospatial/delaunay/data/singapore_train_stations.csv b/content/algorithms/geospatial/delaunay/data/singapore_train_stations.csv
new file mode 100644
index 000000000..8f6a3fa18
--- /dev/null
+++ b/content/algorithms/geospatial/delaunay/data/singapore_train_stations.csv
@@ -0,0 +1,158 @@
+station_name,type,lat,lng
+Jurong East,MRT,1.333207,103.742308
+Bukit Batok,MRT,1.349069,103.749596
+Bukit Gombak,MRT,1.359043,103.751863
+Choa Chu Kang,MRT,1.385417,103.744316
+Yew Tee,MRT,1.397383,103.747523
+Kranji,MRT,1.425302,103.762049
+Marsiling,MRT,1.432579,103.77415
+Woodlands,MRT,1.436984,103.786406
+Admiralty,MRT,1.436984,103.786406
+Sembawang,MRT,1.449133,103.82006
+Yishun,MRT,1.429666,103.835044
+Khatib,MRT,1.417423,103.832995
+Yio Chu Kang,MRT,1.381765,103.844923
+Ang Mo Kio,MRT,1.370025,103.849588
+Bishan,MRT,1.35092,103.848206
+Braddell,MRT,1.34055,103.847098
+Toa Payoh,MRT,1.332405,103.847436
+Novena,MRT,1.320089,103.843405
+Newton,MRT,1.31383,103.838021
+Orchard,MRT,1.304041,103.831792
+Somerset,MRT,1.300508,103.838428
+Dhoby Ghaut,MRT,1.299169,103.845799
+City Hall,MRT,1.293119,103.852089
+Raffles Place,MRT,1.284001,103.85155
+Marina Bay,MRT,1.276481,103.854598
+Marina South Pier,MRT,1.271422,103.863581
+Tuas Link,MRT,1.340371,103.636866
+Tuas West Road,MRT,1.330075,103.639636
+Tuas Crescent,MRT,1.321091,103.649075
+Gul Circle,MRT,1.319809,103.66083
+Joo Koon,MRT,1.327826,103.678318
+Pioneer,MRT,1.337645,103.69742
+Boon Lay,MRT,1.33862,103.705817
+Lakeside,MRT,1.344264,103.720797
+Chinese Garden,MRT,1.342436,103.732582
+Clementi,MRT,1.314925,103.765341
+Dover,MRT,1.311414,103.778596
+Buona Vista,MRT,1.307337,103.790046
+Commonwealth,MRT,1.302439,103.798326
+Queenstown,MRT,1.294867,103.805902
+Redhill,MRT,1.289674,103.816787
+Tiong Bahru,MRT,1.286555,103.826956
+Outram Park,MRT,1.280319,103.839459
+Tanjong Pagar,MRT,1.276385,103.846771
+Bugis,MRT,1.300747,103.855873
+Lavender,MRT,1.307577,103.863155
+Kallang,MRT,1.311532,103.871372
+Aljunied,MRT,1.316474,103.882762
+Paya Lebar,MRT,1.318214,103.893133
+Eunos,MRT,1.319809,103.902888
+Kembangan,MRT,1.320998,103.913433
+Bedok,MRT,1.324043,103.930205
+Tanah Merah,MRT,1.327309,103.946479
+Simei,MRT,1.343237,103.953343
+Tampines,MRT,1.354467,103.943325
+Pasir Ris,MRT,1.373234,103.949343
+Expo,MRT,1.334479,103.961459
+Changi Airport,MRT,1.357622,103.988487
+HarbourFront,MRT,1.265453,103.820514
+Chinatown,MRT,1.284566,103.843626
+Clarke Quay,MRT,1.288949,103.847521
+Little India,MRT,1.306691,103.849396
+Farrer Park,MRT,1.312679,103.854872
+Boon Keng,MRT,1.320091,103.861655
+Potong Pasir,MRT,1.331316,103.868779
+Woodleigh,MRT,1.339202,103.870727
+Serangoon,MRT,1.349862,103.873635
+Kovan,MRT,1.360207,103.885163
+Hougang,MRT,1.371406,103.892533
+Buangkok,MRT,1.382991,103.893347
+Sengkang,MRT,1.391682,103.895475
+Punggol,MRT,1.405191,103.902367
+Bras Basah,MRT,1.296978,103.850715
+Esplanade,MRT,1.293995,103.855396
+Promenade,MRT,1.294063,103.860156
+Nicoll Highway,MRT,1.300292,103.863449
+Stadium,MRT,1.302847,103.875417
+Mountbatten,MRT,1.306106,103.883175
+Dakota,MRT,1.308474,103.888825
+MacPherson,MRT,1.326769,103.889901
+Tai Seng,MRT,1.33594,103.887706
+Bartley,MRT,1.342923,103.87966
+Lorong Chuan,MRT,1.35153,103.864957
+Marymount,MRT,1.349089,103.839116
+Caldecott,MRT,1.337649,103.839627
+Botanic Gardens,MRT,1.322387,103.814905
+Farrer Road,MRT,1.317606,103.807711
+Holland Village,MRT,1.311189,103.796119
+one-north,MRT,1.299854,103.787584
+Kent Ridge,MRT,1.293629,103.784441
+Haw Par Villa,MRT,1.283149,103.781991
+Pasir Panjang,MRT,1.276111,103.791893
+Labrador Park,MRT,1.27218,103.802557
+Telok Blangah,MRT,1.270769,103.809878
+Bayfront,MRT,1.281371,103.858998
+Bukit Panjang,MRT,1.37834,103.762452
+Cashew,MRT,1.369997,103.764569
+Hillview,MRT,1.363185,103.767371
+Beauty World,MRT,1.341607,103.775682
+King Albert Park,MRT,1.335721,103.783203
+Sixth Avenue,MRT,1.331221,103.79718
+Tan Kah Kee,MRT,1.325826,103.807959
+Stevens,MRT,1.320012,103.825964
+Rochor,MRT,1.303601,103.852581
+Downtown,MRT,1.27949,103.852802
+Telok Ayer,MRT,1.282285,103.848584
+Fort Canning,MRT,1.291631,103.844621
+Bencoolen,MRT,1.298477,103.849984
+Jalan Besar,MRT,1.305551,103.855443
+Bendemeer,MRT,1.313674,103.863098
+Geylang Bahru,MRT,1.321479,103.871457
+Mattar,MRT,1.326878,103.883304
+Ubi,MRT,1.330008,103.898911
+Kaki Bukit,MRT,1.335076,103.909057
+Bedok North,MRT,1.335268,103.918054
+Bedok Reservoir,MRT,1.336595,103.93307
+Tampines West,MRT,1.345583,103.938244
+Tampines East,MRT,1.35631,103.955471
+Upper Changi,MRT,1.342218,103.961505
+South View,LRT,1.380299,103.745286
+Keat Hong,LRT,1.378604,103.749058
+Teck Whye,LRT,1.376738,103.753665
+Phoenix,LRT,1.378798,103.758021
+Senja,LRT,1.382852,103.762312
+Jelapang,LRT,1.386703,103.764547
+Segar,LRT,1.387713,103.769599
+Fajar,LRT,1.384502,103.770862
+Bangkit,LRT,1.380281,103.772576
+Pending,LRT,1.376223,103.771277
+Petir,LRT,1.377828,103.76655
+Compassvale,LRT,1.394615,103.900443
+Rumbia,LRT,1.391553,103.905947
+Bakau,LRT,1.38804,103.905412
+Kangkar,LRT,1.383957,103.90216
+Ranggung,LRT,1.384116,103.897386
+Cheng Lim,LRT,1.396332,103.89379
+Farmway,LRT,1.397178,103.889168
+Kupang,LRT,1.398271,103.881283
+Thanggam,LRT,1.397378,103.87561
+Fernvale,LRT,1.392033,103.876256
+Layar,LRT,1.392141,103.880022
+Tongkang,LRT,1.389519,103.885829
+Renjong,LRT,1.386827,103.890541
+Cove,LRT,1.399534,103.905792
+Meridian,LRT,1.397002,103.908884
+Coral Edge,LRT,1.39392,103.912633
+Riviera,LRT,1.39454,103.916056
+Kadaloor,LRT,1.399633,103.916536
+Oasis,LRT,1.402304,103.912736
+Damai,LRT,1.405293,103.908606
+Sam Kee,LRT,1.409808,103.90492
+Teck Lee,LRT,1.412783,103.906565
+Punggol Point,LRT,1.416932,103.90668
+Samudera,LRT,1.415955,103.902185
+Nibong,LRT,1.411865,103.900321
+Sumang,LRT,1.408501,103.898605
+Soo Teck,LRT,1.405436,103.897287
diff --git a/content/algorithms/geospatial/delaunay/delaunay.md b/content/algorithms/geospatial/delaunay/delaunay.md
new file mode 100644
index 000000000..10266165b
--- /dev/null
+++ b/content/algorithms/geospatial/delaunay/delaunay.md
@@ -0,0 +1,276 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.15.0
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+# Delaunay Graph from Geographic Points
+
+Delaunay triangulation for a set P of points in the plane is a triangulation DT(P) such that no point in P is inside the circumcircle of any triangle in DT(P). Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation; they tend to avoid skinny triangles. The triangulation was invented by Boris Delaunay in 1934.
+
+![delaunay_triangulation](images/delaunay_introduction.png "A Delaunay triangulation in the plane with circumcircles shown")
+
+## Properties of Denaulay Triangulation
+
+- Empty Circle Property of Triangles: the circumcircle of $\Delta pqr$ does not contain any other site of $P$.
+- Empty Circle Property of Edges: A pair $(p_i, p_j)$ is an edge of DT if and only if there exists an empty circle passing through $p_i, p_j$.
+- Closest Pair Property: Given a point set $P$, if $p_i, p_j$ are the two closest pair of points, then $(p_i, p_j)$ is an edge of DT.
+- Minimum Spanning Tree: Another nice property of DT is that the minimum spanning tree of the sites is a subgraph of DT.
+- Maximizing the Minimum Angle: In the plane, the Delaunay triangulation maximizes the minimum angle. Compared to any other triangulation of the points, the smallest angle in the Delaunay triangulation is at least as large as the smallest angle in any other. However, the Delaunay triangulation does not necessarily minimize the maximum angle.
+- A circle circumscribing any Delaunay triangle does not contain any other input points in its interior.
+- In the plane, each vertex has on average six surrounding triangles.
+
+## Visual Delaunay definition: Flipping
+
+
+
+![edge_flip](images/edge_flip.png "Edge flipping property of Delaunay Triangulation")
+
+From the above properties an important feature arises: Looking at two triangles **ABD** and **BCD** with the common edge BD (see figures), if the sum of the angles $\alpha$ and  $\gamma$ is less than or equal to $180^{\circ}$, the triangles meet the Delaunay condition.
+
+This is an important property because it allows the use of a flipping technique. If two triangles do not meet the Delaunay condition, switching the common edge BD for the common edge AC produces two triangles that do meet the Delaunay condition.
+
+
+## Algorithms
+
+In two dimensions, one way to detect if point D lies in the circumcircle of A, B, C is to evaluate the determinant:
+
+$$
+\begin{bmatrix}
+A_x & A_y & A_x^2 + A_y^2 & 1 \\
+B_x & B_y & B_x^2 + B_y^2 & 1 \\
+C_x & C_y & C_x^2 + C_y^2 & 1 \\
+D_x & D_y & D_x^2 + D_y^2 & 1 
+\end{bmatrix}
+= 
+\begin{bmatrix}
+A_x - D_x & A_y - D_y & (A_x^2 - D_x^2) + (A_y^2 - D_y^2) \\
+B_x - D_x & B_y - D_y & (B_x^2 - D_x^2) + (B_y^2 - D_y^2) \\
+C_x - D_x & C_y - D_y & (C_x^2 - D_x^2) + (C_y^2 - D_y^2) 
+\end{bmatrix}
+> 0
+$$
+
+Assuming A, B and C to lie counter-clockwise, this is positive if and only if D lies in the circumcircle.
+
+### Flip Algorithms
+
+As mentioned above, if a triangle is non-Delaunay, we can flip one of its edges. This leads to a straightforward algorithm: construct any triangulation of the points, and then flip edges until no triangle is non-Delaunay. Unfortunately, this can take $O(n^2)$ edge flips, and does not extend to three dimensions or higher.
+
+### Incremental
+
+The most straightforward way of efficiently computing the Delaunay triangulation is to repeatedly add one vertex at a time, retriangulating the affected parts of the graph. When a vertex $p$ is added (Fig.a), we split in three the triangle that contains $p$, then we apply the flip algorithm (Fig.b). Done naively, this will take $O(n)$ time: we search through all the triangles to find the one that contains $p$, then we potentially flip away every triangle. Then the overall runtime is $O(n^2)$.
+
+![incremental](images/incremental.png "Incremental Approach in generating Delaunay Triangles")
+
+
+### Divide and Conquer
+
+In this algorithm, one recursively draws a line to split the vertices into two sets. The Delaunay triangulation is computed for each set, and then the two sets are merged along the splitting line. Using some clever tricks, the merge operation can be done in time $O(n)$, so the total running time is $O(nlogn)$.
+
+## Applications
+
+Connecting the centers of the circumcircles in Delaunay triangulation forms the Voronoi diagram. Voronoi diagram and the triangulations are useful in various geometrical applications.
+
+- A point location data structure can be built on top of the Voronoi diagram in order to answer nearest neighbor queries. For example, when one wants to find the nearest hospital.
+- One can also find the largest empty circle amongst a set of points, and in an enclosing polygon. For example, to build a new supermarket as far as possible from all the existing ones.
+- Useful in polymer physics. It can be used to represent free volume of the polymer.
+- In climatology, Voronoi diagrams are used to calculate the rainfall of an area, based on a series of point measurements. In this usage, they are generally referred to as Thiessen polygons.
+- Used to study the growth patterns of forests and forest canopies, and may also be helpful in developing predictive models for forest fires.
+
+## Delaunay Triangulation implementation
+
+Delaunay triangles alongwith the Voronoi diagram can be implemented in python using [PySAL](https://pysal.org/). It provides a rich suite of spatial analysis algorithms.
+
+The `voronoi_frames` method can be used to form the required operation. The following examples will show us how to build delaunay graph (with Voronoi diagrams) from a set of points.
+
+For basic understanding, we will use a set of random points to elaborate the usage of the function and its various parameters.
+
+- Necessary libraries are imported
+
+```{code-cell}
+from libpysal import weights, examples
+from libpysal.cg.voronoi import voronoi_frames, voronoi
+from shapely.geometry import Point
+import matplotlib.pyplot as plt
+import networkx as nx
+import numpy as np
+```
+
+- Declaring the points to be plotted in delaunay graph
+
+```{code-cell}
+points = [
+    (10.2, 5.1),
+    (4.7, 2.2),
+    (5.3, 5.7),
+    (2.7, 5.3),
+    (10.3, 6.4),
+    (4.5, 7),
+    (20, 10),
+    (0, 2),
+]
+regions, vertices = voronoi(points)
+```
+
+- Printing the regions and points of the triangulations
+
+```{code-cell}
+regions
+```
+
+```{code-cell}
+vertices
+```
+
+### Visualizing the Points
+
+The `voronoi_frames` takes in four parameters:
+
+- **points** (array): The originator points
+- **radius** (float): The distance to ‘points at infinity’ used in building voronoi cells. Default is None.
+- **clip** (string): An overloaded option about how to clip the voronoi cells. Default is 'extent'.
+- **tolerance** (float): The percent of map width to use to buffer the extent of the map (default: .01, or 1%).
+
+Next, plotting the aformentioned points in a delaunay graph
+
+```{code-cell}
+region_df, point_df = voronoi_frames(points, clip="extent")
+fig, ax = plt.subplots()
+ax.set_title("Basic example of voronoi_frames", fontsize=12)
+region_df.plot(ax=ax, color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax, color="red");
+```
+
+- Plotting the graph for a radius of 5.0. Here, the radius refers to the distance to ‘points at infinity’ used in building voronoi cells. The default value of radius is `None`.
+
+```{code-cell}
+fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 6))
+
+region_df, point_df = voronoi_frames(points, clip="none", radius=0.0)
+ax[0].set_title("Voronoi Diagram with radius value of 0.0", fontsize=12)
+region_df.plot(ax=ax[0], color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax[0], color="red")
+
+region_df, point_df = voronoi_frames(points, clip="none", radius=100.0)
+ax[1].set_title("Voronoi Diagram with radius value of 100.0", fontsize=12)
+region_df.plot(ax=ax[1], color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax[1], color="red");
+```
+
+Finally, generating graph for different clip options.The clip options available are - 
+- `'none'` / `None`: No clip is applied
+- `'bbox'` / `'extent'` / `'bounding box'`: Clip the voronoi cells to the bounding box of the input points
+- `'chull'` / `'convex hull'`: Clip the voronoi cells to the convex hull of the input points
+- `'ashape'` / `'ahull'`: Clip the voronoi cells to the tightest hull that contains all points
+
+```{code-cell}
+fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(12, 9))
+
+region_df, point_df = voronoi_frames(points, clip="none")
+ax[0, 0].set_title("Voronoi Diagram with no clipping", fontsize=12)
+region_df.plot(ax=ax[0, 0], color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax[0, 0], color="red")
+
+region_df, point_df = voronoi_frames(points, clip="extent")
+ax[0, 1].set_title("Voronoi Diagram with Bounding Box clipping", fontsize=12)
+region_df.plot(ax=ax[0, 1], color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax[0, 1], color="red")
+
+region_df, point_df = voronoi_frames(points, 100.0, clip="chull")
+ax[1, 0].set_title("Voronoi Diagram with Convex Hull clipping", fontsize=12)
+region_df.plot(ax=ax[1, 0], color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax[1, 0], color="red")
+
+region_df, point_df = voronoi_frames(points, 100.0, clip="ahull")
+ax[1, 1].set_title("Voronoi Diagram with tightest hull clipping", fontsize=12)
+region_df.plot(ax=ax[1, 1], color="blue", edgecolor="black", alpha=0.3)
+point_df.plot(ax=ax[1, 1], color="red");
+```
+
+### Voronoi Frames with Real World Dataset
+
+This example shows how to build a delaunay graph (plus its dual, the set of Voronoi polygons) from a set of points. For this, we will use the dataset of train stations in Singapore. The methods shown here can also work directly with polygonal data using their centroids as representative points.
+
+- First, including the necessary libraries
+
+```{code-cell}
+from libpysal import weights, examples
+from libpysal.cg import voronoi_frames
+from contextily import add_basemap
+import matplotlib.pyplot as plt
+import networkx as nx
+import pandas as pd
+import numpy as np
+```
+
+- Reading the dataset from a CSV file using Pandas. The dataset contains the information about each train stations in Singapore alongwith its latitude and longitudes. Next, the coordinates are converted into float data type from string
+
+```{code-cell}
+data = pd.read_csv("data/singapore_train_stations.csv")
+
+data.lng = list(map(float, data.lng))
+data.lat = list(map(float, data.lat))
+```
+
+- In order for networkx to plot the nodes of our graph correctly, we need to construct the array of coordinates for each point in our dataset. To get this as a numpy array, we extract the x and y coordinates from the geometry column.
+
+```{code-cell}
+coordinates = np.column_stack((data.lng, data.lat))
+```
+
+- While we could simply present the Delaunay graph directly, it is useful to visualize the Delaunay graph alongside the Voronoi diagram. This is because the two are intrinsically linked: the adjacency graph of the Voronoi diagram is the Delaunay graph for the set of generator points! Put simply, this means we can build the Voronoi diagram (relying on scipy.spatial for the underlying computations), and then convert these polygons quickly into the Delaunay graph. Be careful, though; our algorithm, by default, will clip the voronoi diagram to the bounding box of the point pattern. This is controlled by the "clip" argument.
+
+```{code-cell}
+cells, generators = voronoi_frames(coordinates, clip="chull")
+```
+
+- With the voronoi polygons, we can construct the adjacency graph between them using "Rook" contiguity. This represents voronoi cells as being adjacent if they share an edge/face. This is an analogue to the "von Neuman" neighborhood, or the 4 cardinal neighbors in a regular grid. The name comes from the directions a Rook piece can move on a chessboard.
+
+```{code-cell}
+delaunay = weights.Rook.from_dataframe(cells)
+```
+
+- Once the graph is built, we can convert the graphs to networkx objects using the relevant method.
+
+```{code-cell}
+delaunay_graph = delaunay.to_networkx()
+```
+
+- To plot with networkx, we need to merge the nodes back to their positions in order to plot in networkx
+
+```{code-cell}
+positions = dict(zip(delaunay_graph.nodes, coordinates))
+```
+
+- Now, we can plot with a nice basemap using Contextily
+
+```{code-cell}
+figure, ax = plt.subplots(figsize=(15, 9))
+cells.plot(ax=ax, facecolor="lightblue", alpha=0.50, edgecolor="cornsilk", linewidth=2)
+add_basemap(ax, crs="EPSG:4326")
+ax.axis("off")
+nx.draw(
+    delaunay_graph,
+    positions,
+    ax=ax,
+    node_size=2,
+    node_color="k",
+    edge_color="k",
+    alpha=0.8,
+)
+plt.show()
+```
+
+## References
+
+- [Delaunay Triangulation](http://wiki.gis.com/wiki/index.php/Delaunay_triangulation)
+- [libpysal.cg.voronoi_frames](https://pysal.org/libpysal/generated/libpysal.cg.voronoi_frames.html#libpysal.cg.voronoi_frames)
+- [Incremental approach](https://www.cs.umd.edu/class/spring2020/cmsc754/Lects/lect13-delaun-alg.pdf)
diff --git a/content/algorithms/geospatial/delaunay/images/delaunay_introduction.png b/content/algorithms/geospatial/delaunay/images/delaunay_introduction.png
new file mode 100644
index 000000000..ac0111028
Binary files /dev/null and b/content/algorithms/geospatial/delaunay/images/delaunay_introduction.png differ
diff --git a/content/algorithms/geospatial/delaunay/images/edge_flip.png b/content/algorithms/geospatial/delaunay/images/edge_flip.png
new file mode 100644
index 000000000..6416a388a
Binary files /dev/null and b/content/algorithms/geospatial/delaunay/images/edge_flip.png differ
diff --git a/content/algorithms/geospatial/delaunay/images/incremental.png b/content/algorithms/geospatial/delaunay/images/incremental.png
new file mode 100644
index 000000000..8e1e6d76c
Binary files /dev/null and b/content/algorithms/geospatial/delaunay/images/incremental.png differ
diff --git a/content/algorithms/geospatial/index.md b/content/algorithms/geospatial/index.md
new file mode 100644
index 000000000..005decd74
--- /dev/null
+++ b/content/algorithms/geospatial/index.md
@@ -0,0 +1,11 @@
+# Geospatial Information System
+
+A closer look at techniques and libraries recommended by Networkx when working with geospatial data (including reading and writing shapefiles).
+
+```{toctree}
+---
+maxdepth: 1
+---
+delaunay/delaunay
+
+```
diff --git a/requirements.txt b/requirements.txt
index 6fa8a0337..b68a132cd 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -20,3 +20,8 @@ pyupgrade
 # Notebook requirements
 pygraphviz
 
+libpysal
+contextily
+geopandas
+shapely
+rtree
diff --git a/site/index.md b/site/index.md
index 319059d54..1bad56b99 100644
--- a/site/index.md
+++ b/site/index.md
@@ -36,6 +36,7 @@ maxdepth: 1
 ---
 
 content/algorithms/index
+content/algorithms/geospatial/index
 content/generators/index
 content/exploratory_notebooks/facebook_notebook
 ```