diff --git a/growbikenet/functions.py b/growbikenet/functions.py
index 9cef472..4422f5a 100644
--- a/growbikenet/functions.py
+++ b/growbikenet/functions.py
@@ -28,22 +28,25 @@ def intersects_properly(geom1, geom2):
     return geom1.intersects(geom2) and not geom1.touches(geom2)
 
 
-def prepare_network(city_name, proj_crs, network_type='all_public', custom_filter=None, retain_all=True):
+def prepare_network(city_name, proj_crs, network_type='all_public', custom_filter=None, retain_all=True, city_boundary_file=None):
     """Download and prepare a street network from OSM via OSMnx
     Downloads a network with a given network_type and custom_filter using ox.graph_from_place.
     Then, stores the undirected OSM data in gdfs and projects using proj_crs.
     Parameters
     ----------
     city_name : str
-        Name of the city that the analysis should be performed on.
-    proj_crs : str, default '3857'
-        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+        Name of the city that the analysis should be performed on. Overruled (for data fetching) if city_boundary_file is set.
+    proj_crs : str
+        Coordinate reference system that is used to project osm data.
     network_type : {“all”, “all_public”, “bike”, “drive”, “drive_service”, “walk”} 
         What type of street network to retrieve if custom_filter is None.
     custom_filter : (str | list[str] | None)
         A custom ways filter to be used instead of the network_type presets
     retain_all : bool, default True
         If True, return the entire graph even if it is not connected, useful for disconnected bicycle networks. If False, retain only the largest weakly connected component, useful for road networks.
+    city_boundary_file : (str | None), default None
+        If not set to None, the study area will be selected from the (Multi)Polygon provided in the city_boundary_file shape file. For example, "copenhagen.shp".
+
     Returns
     -------
     nodes : geopandas.geodataframe.GeoDataFrame
@@ -53,25 +56,34 @@ def prepare_network(city_name, proj_crs, network_type='all_public', custom_filte
     g_undir : networkx.classes.multigraph.MultiGraph
         Extracted networkX graph, undirected
     """
+
     # Fetch street network data from osmnx
-    g = ox.graph_from_place(
-    city_name, network_type=network_type, custom_filter=custom_filter, retain_all=retain_all
-    )
+    if city_boundary_file is None:
+        g = ox.graph_from_place(
+        city_name, network_type=network_type, custom_filter=custom_filter, retain_all=retain_all
+        )
+    else:
+        shp = gpd.read_file(city_boundary_file)
+        city_boundary_polygon = shp.iloc[0].geometry
+        g = ox.graph_from_polygon(
+        city_boundary_polygon, network_type=network_type, custom_filter=custom_filter, retain_all=retain_all
+        )
+
     g_undir = g.to_undirected().copy() # convert to undirected (dropping OSMnx keys!)
 
     # Export osmnx data to gdfs
     nodes, edges = nx_to_nodes_edges(g_undir, proj_crs)
     return nodes, edges, g_undir
 
-def nx_to_nodes_edges(G, proj_crs='3857'):
+def nx_to_nodes_edges(G, proj_crs):
     """Get nodes and projected edges from networkX graph
     
     Parameters
     ----------
     G : networkx.classes.multigraph.MultiGraph
         networkX graph, undirected
-    proj_crs : str, default '3857'
-        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+    proj_crs : str
+        Coordinate reference system that is used to project osm data.
         
     Returns
     -------
@@ -165,7 +177,7 @@ def get_existing_network_seed_points(nodes_exnw, existing_network_spacing):
 
     return seed_points_exnw
     
-def update_with_existing_bike_network(city_name, proj_crs, g_undir):
+def update_with_existing_bike_network(city_name, proj_crs, g_undir, city_boundary_file=None):
     """Update street network with existing bike network
 
     Downloads a network of protected bike infrastructure from OSM (retaining all connected components) and merges it to a given street network graph g_undir.
@@ -173,11 +185,13 @@ def update_with_existing_bike_network(city_name, proj_crs, g_undir):
     Parameters
     ----------
     city_name : str
-        Name of the city that the analysis should be performed on.
+        Name of the city that the analysis should be performed on. Overruled (for data fetching) if city_boundary_file is set.
     proj_crs : str
-        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+        Coordinate reference system that is used to project osm data.
     g_undir : networkx.classes.multigraph.MultiGraph
         Street network networkX graph, undirected
+    city_boundary_file : (str | None), default None
+        If not set to None, the study area will be selected from the (Multi)Polygon provided in the city_boundary_file shape file. For example, "copenhagen.shp".
 
     Returns
     -------
@@ -205,7 +219,7 @@ def update_with_existing_bike_network(city_name, proj_crs, g_undir):
             ox.settings.useful_tags_way.extend(custom_tag)
     # Fetch protected bike network data from osmnx
     # Due to retain_all=True, this fetches all the connected components
-    nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf, retain_all=True)
+    nodes_exnw, edges_exnw, g_undir_exnw = prepare_network(city_name, proj_crs, custom_filter=cf, retain_all=True, city_boundary_file=city_boundary_file)
     g_undir = nx.compose(g_undir_exnw, g_undir) # Merge to be sure we have everything from both
 
     # Now we could have some leftover bike infra that is disconnected from the street network and thus not routable.
@@ -237,7 +251,7 @@ def update_seed_points_with_existing_bike_network(seed_points_snapped, nodes_exn
     existing_network_spacing : int
         Positive integer denoting spacing between seed points, in meters, only on the existing bicycle network.
     proj_crs : str
-        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator).
+        Coordinate reference system that is used to project osm data.
 
     Returns
     -------
diff --git a/growbikenet/growbikenet.py b/growbikenet/growbikenet.py
index b64f088..3977864 100644
--- a/growbikenet/growbikenet.py
+++ b/growbikenet/growbikenet.py
@@ -38,6 +38,7 @@ def growbikenet(
     export_plots=False,
     export_video=False,
     allow_edge_overlaps=False,
+    city_boundary_file=None,
 ):
     """Creates a list of edges ordered by a specified ranking method.
 
@@ -46,9 +47,9 @@ def growbikenet(
     Parameters
     ----------
     city_name : str
-        Name of the city that the analysis should be performed on
+        Name of the city that the analysis should be performed on. Overruled (for data fetching) if city_boundary_file is set.
     proj_crs : str, default '3857'
-        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator)
+        Coordinate reference system that is used to project osm data. Default is '3857' (WGS 84 / Pseudo-Mercator). If this web mercator projection is not needed, then for Europe '3035' (LAEA) and globally '54035' (Equal Earth) is better.
     ranking : str, default 'betweenness_centrality'
         Method used to rank edges. Must be 'betweenness_centrality' (default), 'closeness_centrality', or 'random'.
     seed_point_type : str, optional, default 'grid'
@@ -62,7 +63,7 @@ def growbikenet(
         Spacing between seed points, in meters, only on the existing bicycle network. If not set to a positive integer, the existing network is ignored.
     export_data : bool, optional, default True
         If set to True, data will be saved to a file. The filename is [slug]-[ranking]-[seed_point_type].gpkg, where slug is a string id made out of city_name
-    export_data_slug : stri, optional, default None
+    export_data_slug : str, optional, default None
         If not set to None, the city_name will be slugified and used as the slug in the filename of the data export
     export_file_format : str, optional, default "geojson"
         File format for the data export, relevant if export_data set to True. Default "geojson", also possible "gpkg". If exporting as geojson, generates extra files for seed points and city boundary. If exporting as gkpg, these are added all in one file as extra layers.
@@ -72,6 +73,8 @@ def growbikenet(
         If set to True, video will be saved to a file (only possible if export_plots is set to True)
     allow_edge_overlaps : bool, default False
         If set to False, removes edge overlaps in consecutive growth stages. In this case, growth stages that do not add anything new are deleted.
+    city_boundary_file : (str | None), default None
+        If not set to None, the study area will be selected from the (Multi)Polygon provided in the city_boundary_file shape file, ideally in unprojected latitude-longitude degrees (EPSG:4326), but EPSG:3857 also works. For example, "./tests/test_data/copenhagen.shp".
 
     Returns
     -------
@@ -84,7 +87,8 @@ def growbikenet(
     .. [2] P. Folco, L. Gauvin, M. Tizzoni, M. Szell, "Data-driven micromobility network planning for demand and safety", Environment and planning B: Urban analytics and city science 50(8), 2087-2102 (2023)
 
     """
-    # check if user input is valid
+    
+    # Check if user input is valid
     if type(city_name) is not str:
         raise TypeError("city_name must be a string")
     if type(proj_crs) is not str:
@@ -129,18 +133,26 @@ def growbikenet(
         raise TypeError("export_plots must be a boolean")
     if type(export_video) is not bool:
         raise TypeError("export_video must be a boolean")
+    if city_boundary_file is not None and type(city_boundary_file) is not str:
+        raise TypeError("city_boundary_file must be None or a string")
+    if type(city_boundary_file) is str and not os.path.isfile(city_boundary_file):
+        raise FileNotFoundError("city_boundary_file not found")
 
     np.random.seed(42)  # Set random number generator seed for reproducibility
 
     ### Download and preprocess data from OSM
-    print("Downloading OSM data..")
+    if city_boundary_file:
+        import_string = " from city boundary provided in "+city_boundary_file
+    else:
+        import_string = ""
+    print("Downloading OSM data"+import_string+"..")
 
     # Fetch street network data from osmnx
     # Due to retain_all=False, this fetches the largest connected component
-    nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public', retain_all=False)
+    nodes, edges, g_undir = prepare_network(city_name, proj_crs, network_type='all_public', retain_all=False, city_boundary_file=city_boundary_file)
 
     if existing_network_spacing: # TO DO: Check for empty bike infra!
-        nodes, edges, g_undir, nodes_exnw, edges_exnw = update_with_existing_bike_network(city_name, proj_crs, g_undir)
+        nodes, edges, g_undir, nodes_exnw, edges_exnw = update_with_existing_bike_network(city_name, proj_crs, g_undir, city_boundary_file=city_boundary_file)
 
     ### Create seed points
     print("Creating " + seed_point_type + " seed points..")
@@ -270,7 +282,11 @@ def growbikenet(
         ### save data
         print("Saving data..")
         seed_points_snapped.drop(["osmid"], axis=1, inplace=True)
-        city_boundary = ox.geocoder.geocode_to_gdf(city_name)
+        if city_boundary_file:
+            shp = gpd.read_file(city_boundary_file)
+            city_boundary = shp.iloc[[0]] # This needs to stay a gdf
+        else:
+            city_boundary = ox.geocoder.geocode_to_gdf(city_name)
         city_boundary.to_crs(epsg=proj_crs, inplace=True)
         # We have meter precision, so rounding to integers is fine. Better would be to 
         # change dtypes to int, but this does not seem possible without manual looping.
diff --git a/tests/test_data/copenhagen.dbf b/tests/test_data/copenhagen.dbf
new file mode 100644
index 0000000..101a0ec
Binary files /dev/null and b/tests/test_data/copenhagen.dbf differ
diff --git a/tests/test_data/copenhagen.prj b/tests/test_data/copenhagen.prj
new file mode 100644
index 0000000..a30c00a
--- /dev/null
+++ b/tests/test_data/copenhagen.prj
@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
\ No newline at end of file
diff --git a/tests/test_data/copenhagen.qpj b/tests/test_data/copenhagen.qpj
new file mode 100644
index 0000000..5fbc831
--- /dev/null
+++ b/tests/test_data/copenhagen.qpj
@@ -0,0 +1 @@
+GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
diff --git a/tests/test_data/copenhagen.shp b/tests/test_data/copenhagen.shp
new file mode 100644
index 0000000..81dc6e6
Binary files /dev/null and b/tests/test_data/copenhagen.shp differ
diff --git a/tests/test_data/copenhagen.shx b/tests/test_data/copenhagen.shx
new file mode 100644
index 0000000..66e4f92
Binary files /dev/null and b/tests/test_data/copenhagen.shx differ