diff --git a/examples/duplicate_points.png b/examples/duplicate_points.png
index ef5cdce..45c9471 100644
Binary files a/examples/duplicate_points.png and b/examples/duplicate_points.png differ
diff --git a/examples/random_points_across_italy.png b/examples/random_points_across_italy.png
index b01369b..ec08af8 100644
Binary files a/examples/random_points_across_italy.png and b/examples/random_points_across_italy.png differ
diff --git a/examples/random_points_across_italy_per_geom_false.png b/examples/random_points_across_italy_per_geom_false.png
index d44054e..bc63b2e 100644
Binary files a/examples/random_points_across_italy_per_geom_false.png and b/examples/random_points_across_italy_per_geom_false.png differ
diff --git a/examples/random_points_and_area.png b/examples/random_points_and_area.png
index 7321b49..4077819 100644
Binary files a/examples/random_points_and_area.png and b/examples/random_points_and_area.png differ
diff --git a/examples/random_points_brandenburg.png b/examples/random_points_brandenburg.png
index b916499..416f203 100644
Binary files a/examples/random_points_brandenburg.png and b/examples/random_points_brandenburg.png differ
diff --git a/examples/using_geopandas.png b/examples/using_geopandas.png
index 948f250..3842546 100644
Binary files a/examples/using_geopandas.png and b/examples/using_geopandas.png differ
diff --git a/examples/using_geopandas.py b/examples/using_geopandas.py
index 609b73f..67a3dad 100644
--- a/examples/using_geopandas.py
+++ b/examples/using_geopandas.py
@@ -9,7 +9,7 @@
 
 import matplotlib.pyplot as plt
 import geopandas as gpd
-from shapely.ops import cascaded_union
+from shapely.ops import unary_union
 
 from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
 from geovoronoi import voronoi_regions_from_coords, points_to_coords
@@ -34,7 +34,7 @@
 cities = cities.to_crs(south_am.crs)   # convert city coordinates to same CRS!
 
 # create the bounding shape as union of all South American countries' shapes
-south_am_shape = cascaded_union(south_am.geometry)
+south_am_shape = unary_union(south_am.geometry)
 south_am_cities = cities[cities.geometry.within(south_am_shape)]   # reduce to cities in South America
 
 #%%
diff --git a/geovoronoi/_voronoi.py b/geovoronoi/_voronoi.py
index 70bda96..83ea26e 100644
--- a/geovoronoi/_voronoi.py
+++ b/geovoronoi/_voronoi.py
@@ -10,11 +10,15 @@
 from collections import defaultdict
 
 import numpy as np
-from scipy.spatial import Voronoi
-from scipy.spatial.qhull import QhullError
-from shapely.geometry import box, LineString, asPoint, MultiPoint, Polygon, MultiPolygon
+
+try:
+    from scipy.spatial import Voronoi, QhullError
+except ImportError:
+    from scipy.spatial.qhull import QhullError    # for older versions of scipy
+
+from shapely.geometry import box, LineString, Point, MultiPoint, Polygon, MultiPolygon
 from shapely.errors import TopologicalError
-from shapely.ops import cascaded_union
+from shapely.ops import unary_union
 
 from ._geom import line_segment_intersection
 
@@ -32,7 +36,7 @@ def coords_to_points(coords):
     :param coords: NumPy array of shape (N,2) with N coordinates in 2D space
     :return: list of length N with Shapely Point objects
     """
-    return list(map(asPoint, coords))
+    return list(map(Point, coords))
 
 
 def points_to_coords(pts):
@@ -276,7 +280,7 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
     geom_bb = box(*geom.bounds).buffer(bounds_buf)
 
     # center of all points
-    center = np.array(MultiPoint(vor.points).convex_hull.centroid)
+    center = np.array(MultiPoint(vor.points).convex_hull.centroid.coords)
 
     # ridge vertice's coordinates
     ridge_vert = np.asarray(vor.ridge_vertices)
@@ -361,7 +365,10 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
                                 isects.append(isect)
 
                         if len(isects) == 0:
-                            raise RuntimeError('ridge line must intersect with surrounding geometry from `geom`')
+                            raise RuntimeError('ridge line must intersect with surrounding geometry from `geom`; '
+                                               'this error often arises when there are points outside of the '
+                                               'surrounding geometries; first check if all your points are inside '
+                                               'the surrounding geometries')
                         elif len(isects) == 1:   # one intersection
                             far_pt = isects[0]
                         else:                    # multiple intersections - take closest intersection
@@ -440,14 +447,14 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
             if len(inner_regions) == 1:
                 inner_regions_poly = region_polys[next(iter(inner_regions))]
             else:
-                inner_regions_poly = cascaded_union([region_polys[i_reg] for i_reg in inner_regions])
+                inner_regions_poly = unary_union([region_polys[i_reg] for i_reg in inner_regions])
 
         # generate polygon from all other regions' polygons other than the current region `i_reg`
         other_regions_polys = [region_polys[i_other] for i_other in border_regions if i_reg != i_other]
         if inner_regions_poly:
             other_regions_polys.append(inner_regions_poly)
 
-        union_other_regions = cascaded_union(other_regions_polys)
+        union_other_regions = unary_union(other_regions_polys)
 
         # generate difference between geom and other regions' polygons -- what's left is the current region's area
         # plus any so far uncovered area
@@ -489,7 +496,7 @@ def region_polygons_from_voronoi(vor, geom, return_point_assignments=False,
         # add new areas as union
         if add:
             old_reg_area = region_polys[i_reg].area
-            new_reg = cascaded_union([p] + add)
+            new_reg = unary_union([p] + add)
             area_diff = new_reg.area - old_reg_area
 
             region_polys[i_reg] = new_reg
diff --git a/setup.py b/setup.py
index 0f8733c..0fc1e7c 100644
--- a/setup.py
+++ b/setup.py
@@ -8,17 +8,17 @@
 GITHUB_URL = 'https://github.com/WZBSocialScienceCenter/geovoronoi'
 
 __title__ = 'geovoronoi'
-__version__ = '0.3.0'
+__version__ = '0.3.1.dev'
 __author__ = 'Markus Konrad'
 __license__ = 'Apache License 2.0'
 
 here = os.path.abspath(os.path.dirname(__file__))
 
-DEPS_BASE = ['numpy>=1.19.0,<2', 'scipy>=1.5.0,<1.7', 'shapely>=1.7.0,<1.8']
+DEPS_BASE = ['numpy>=1.19.0', 'scipy>=1.5.0', 'shapely>=1.7.0']
 
 DEPS_EXTRA = {
-    'plotting': ['matplotlib>=3.3.0,<3.4', 'geopandas>=0.8.0,<0.9', 'descartes>=1.1.0,<1.2'],
-    'test': ['pytest>=6.2.0,<6.3', 'pytest-mpl>=0.12,<0.13', 'hypothesis>=6.0.0,<6.1', 'tox>=3.21.0,<3.22'],
+    'plotting': ['matplotlib>=3.3.0', 'geopandas>=0.8.0', 'descartes>=1.1.0'],
+    'test': ['pytest>=6.2.0', 'pytest-mpl>=0.12', 'hypothesis>=6.0.0', 'tox>=3.21.0'],
     'develop': ['ipython>=7.19.0', 'twine>=3.3.0'],
 }
 
@@ -40,7 +40,7 @@
     url=GITHUB_URL,
     project_urls={
         'Source': GITHUB_URL,
-        'Tracker': 'https://github.com/WZBSocialScienceCenter/geovoronoi' + '/issues',
+        'Tracker': GITHUB_URL + '/issues',
     },
     author='Markus Konrad',
     author_email='markus.konrad@wzb.eu',
@@ -60,6 +60,7 @@
         'Programming Language :: Python :: 3.7',
         'Programming Language :: Python :: 3.8',
         'Programming Language :: Python :: 3.9',
+        'Programming Language :: Python :: 3.10',
 
         'Topic :: Scientific/Engineering :: Information Analysis',
         'Topic :: Software Development :: Libraries :: Python Modules',
diff --git a/tests/test_main.py b/tests/test_main.py
index 2c66b1d..7b10191 100644
--- a/tests/test_main.py
+++ b/tests/test_main.py
@@ -6,8 +6,8 @@
 
 import numpy as np
 import geopandas as gpd
-from shapely.geometry import Polygon, MultiPolygon, asPoint
-from shapely.ops import cascaded_union
+from shapely.geometry import Polygon, MultiPolygon, Point
+from shapely.ops import unary_union
 import pytest
 from hypothesis import given, settings
 import hypothesis.strategies as st
@@ -263,7 +263,7 @@ def test_voronoi_geopandas_with_plot():
     cities = cities.to_crs(south_am.crs)  # convert city coordinates to same CRS!
 
     # create the bounding shape as union of all South American countries' shapes
-    south_am_shape = cascaded_union(south_am.geometry)
+    south_am_shape = unary_union(south_am.geometry)
     south_am_cities = cities[cities.geometry.within(south_am_shape)]  # reduce to cities in South America
 
     # convert the pandas Series of Point objects to NumPy array of coordinates
@@ -407,7 +407,7 @@ def _check_region_polys(region_polys, region_pts, coords, expected_sum_area,
                 polybuf = poly
         else:
             polybuf = poly
-        assert all([polybuf.contains(asPoint(coords[i_pt])) for i_pt in pt_indices])
+        assert all([polybuf.contains(Point(coords[i_pt])) for i_pt in pt_indices])
 
         sum_area += poly.area
 
diff --git a/tox.ini b/tox.ini
index ef0ef24..a1712b1 100644
--- a/tox.ini
+++ b/tox.ini
@@ -4,11 +4,10 @@
 # and then run "tox" from this directory.
 
 [tox]
-envlist = py36, py37, py38, py39
+envlist = py36, py37, py38, py39, py310
 
 [testenv]
 deps = .[test]
 extras = plotting
 commands =
     pytest
-    pytest --mpl --mpl-baseline-path=tests/baseline