diff --git a/aigis/annotate/utils.py b/aigis/annotate/utils.py
index 2fb4b17..ac3d17a 100644
--- a/aigis/annotate/utils.py
+++ b/aigis/annotate/utils.py
@@ -6,6 +6,75 @@
 from pprint import pprint
 from matplotlib import pylab as plt
 import pandas as pd
+from shapely.geometry import Polygon, box
+from pyproj import Proj, transform
+
+
+def create_grid_geojson(bbox, tile_size):
+    min_lon, min_lat, max_lon, max_lat = bbox
+
+    # Transform the coordinates to EPSG:3857 (Web Mercator) for easier calculations
+    in_proj = Proj(init='epsg:4326')
+    out_proj = Proj(init='epsg:3857')
+    min_lon, min_lat = transform(in_proj, out_proj, min_lon, min_lat)
+    max_lon, max_lat = transform(in_proj, out_proj, max_lon, max_lat)
+
+    # Calculate the number of tiles in x and y directions
+    num_tiles_x = int((max_lon - min_lon) / tile_size)
+    num_tiles_y = int((max_lat - min_lat) / tile_size)
+
+    features = []
+
+    for i in range(num_tiles_x):
+        for j in range(num_tiles_y):
+            # Calculate the coordinates of the current tile
+            tile_min_lon = min_lon + i * tile_size
+            tile_max_lon = min_lon + (i + 1) * tile_size
+            tile_min_lat = min_lat + j * tile_size
+            tile_max_lat = min_lat + (j + 1) * tile_size
+
+            # Convert the coordinates back to EPSG:4326
+            tile_min_lon, tile_min_lat = transform(out_proj, in_proj, tile_min_lon, tile_min_lat)
+            tile_max_lon, tile_max_lat = transform(out_proj, in_proj, tile_max_lon, tile_max_lat)
+
+            # Create a polygon for the current tile
+            tile_polygon = box(tile_min_lon, tile_min_lat, tile_max_lon, tile_max_lat)
+
+            # Create a GeoJSON feature for the current tile
+            feature = {
+                "type": "Feature",
+                "properties": {
+                    "id": i * num_tiles_y + j,
+                    "left": tile_min_lon,
+                    "top": tile_max_lat,
+                    "right": tile_max_lon,
+                    "bottom": tile_min_lat,
+                    "row_index": i,
+                    "col_index": j
+                },
+                "geometry": {
+                    "type": "Polygon",
+                    "coordinates": [list(tile_polygon.exterior.coords)]
+                }
+            }
+
+            features.append(feature)
+
+    # Create a GeoJSON feature collection
+    feature_collection = {
+        "type": "FeatureCollection",
+        "name": "GSU_grid_1",
+        "crs": {
+            "type": "name",
+            "properties": {
+                "name": "urn:ogc:def:crs:EPSG::3857"
+            }
+        },
+        "features": features
+    }
+
+    return json.dumps(feature_collection)
+
 
 def show_mask(
     image, mask, alpha=0.1, cmap="viridis", edges=True, edge_colour="green", output=None
diff --git a/scripts/coco2geojson.py b/scripts/coco2geojson.py
index 773e5a1..bbab365 100644
--- a/scripts/coco2geojson.py
+++ b/scripts/coco2geojson.py
@@ -509,5 +509,5 @@ def main(args=None):
     else:
         log.warning("Multiple input tiles detected. The code will handle overlapping and adjacent tiles by default.")
 
-    if __name__ == "__main__":
-        main()
+if __name__ == "__main__":
+    main()
\ No newline at end of file