diff --git a/README.md b/README.md
index 5744d5623..e647ef362 100644
--- a/README.md
+++ b/README.md
@@ -469,6 +469,8 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
  * `-al` or `--drop-lines`: Let "dot" dropping at lower zooms apply to lines too
  * `-ap` or `--drop-polygons`: Let "dot" dropping at lower zooms apply to polygons too
  * `-K` _distance_ or `--cluster-distance=`_distance_: Cluster points (as with `--cluster-densest-as-needed`, but without the experimental discovery process) that are approximately within _distance_ of each other. The units are tile coordinates within a nominally 256-pixel tile, so the maximum value of 255 allows only one feature per tile. Values around 10 are probably appropriate for typical marker sizes. See `--cluster-densest-as-needed` below for behavior.
+ * `-k` _zoom_ or `--cluster-maxzoom=`_zoom_: Max zoom on which to cluster points if clustering is enabled.
+ * `-kg` or `--cluster-maxzoom=g`: Set `--cluster-maxzoom=` to `maxzoom - 1` so that all features are visible at the maximum zoom level.
 
 ### Dropping a fraction of features to keep under tile size limits
 
diff --git a/main.cpp b/main.cpp
index ae41d8254..3da0aa378 100644
--- a/main.cpp
+++ b/main.cpp
@@ -75,6 +75,7 @@ double simplification = 1;
 size_t max_tile_size = 500000;
 size_t max_tile_features = 200000;
 int cluster_distance = 0;
+int cluster_maxzoom = MAX_ZOOM;
 long justx = -1, justy = -1;
 std::string attribute_for_id = "";
 
@@ -1135,7 +1136,7 @@ void choose_first_zoom(long long *file_bbox, std::vector<struct reader> &readers
 	}
 }
 
-int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzoom, int basezoom, double basezoom_marker_width, sqlite3 *outdb, const char *outdir, std::set<std::string> *exclude, std::set<std::string> *include, int exclude_all, json_object *filter, double droprate, int buffer, const char *tmpdir, double gamma, int read_parallel, int forcetable, const char *attribution, bool uses_gamma, long long *file_bbox, const char *prefilter, const char *postfilter, const char *description, bool guess_maxzoom, std::map<std::string, int> const *attribute_types, const char *pgm, std::map<std::string, attribute_op> const *attribute_accum, std::map<std::string, std::string> const &attribute_descriptions, std::string const &commandline) {
+int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzoom, int basezoom, double basezoom_marker_width, sqlite3 *outdb, const char *outdir, std::set<std::string> *exclude, std::set<std::string> *include, int exclude_all, json_object *filter, double droprate, int buffer, const char *tmpdir, double gamma, int read_parallel, int forcetable, const char *attribution, bool uses_gamma, long long *file_bbox, const char *prefilter, const char *postfilter, const char *description, bool guess_maxzoom, bool guess_cluster_maxzoom, std::map<std::string, int> const *attribute_types, const char *pgm, std::map<std::string, attribute_op> const *attribute_accum, std::map<std::string, std::string> const &attribute_descriptions, std::string const &commandline) {
 	int ret = EXIT_SUCCESS;
 
 	std::vector<struct reader> readers;
@@ -1985,7 +1986,7 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 			}
 
 			bool changed = false;
-			while (maxzoom < 32 - full_detail && maxzoom < 33 - low_detail && cluster_distance > 0) {
+			while (maxzoom < 32 - full_detail && maxzoom < 33 - low_detail && maxzoom < cluster_maxzoom && cluster_distance > 0) {
 				unsigned long long zoom_mingap = ((1LL << (32 - maxzoom)) / 256 * cluster_distance) * ((1LL << (32 - maxzoom)) / 256 * cluster_distance);
 				if (avg > zoom_mingap) {
 					break;
@@ -1995,7 +1996,7 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 				changed = true;
 			}
 			if (changed) {
-				printf("Choosing a maxzoom of -z%d to keep most features distinct with cluster distance %d\n", maxzoom, cluster_distance);
+				printf("Choosing a maxzoom of -z%d to keep most features distinct with cluster distance %d and cluster maxzoom %d\n", maxzoom, cluster_distance, cluster_maxzoom);
 			}
 		}
 
@@ -2033,6 +2034,11 @@ int read_input(std::vector<source> &sources, char *fname, int maxzoom, int minzo
 		}
 	}
 
+	if (cluster_maxzoom >= maxzoom && guess_cluster_maxzoom) {
+		cluster_maxzoom = maxzoom - 1;
+		fprintf(stderr, "Choosing a cluster maxzoom of -k%d to make all features visible at maximum zoom %d\n", cluster_maxzoom, maxzoom);
+	}
+
 	if (basezoom < 0 || droprate < 0) {
 		struct tile {
 			unsigned x;
@@ -2478,6 +2484,7 @@ int main(int argc, char **argv) {
 	const char *prefilter = NULL;
 	const char *postfilter = NULL;
 	bool guess_maxzoom = false;
+	bool guess_cluster_maxzoom = false;
 
 	std::set<std::string> exclude, include;
 	std::map<std::string, int> attribute_types;
@@ -2549,6 +2556,7 @@ int main(int argc, char **argv) {
 		{"drop-lines", no_argument, &additional[A_LINE_DROP], 1},
 		{"drop-polygons", no_argument, &additional[A_POLYGON_DROP], 1},
 		{"cluster-distance", required_argument, 0, 'K'},
+		{"cluster-maxzoom", required_argument, 0, 'k'},
 
 		{"Dropping or merging a fraction of features to keep under tile size limits", 0, 0, 0},
 		{"drop-densest-as-needed", no_argument, &additional[A_DROP_DENSEST_AS_NEEDED], 1},
@@ -2795,6 +2803,15 @@ int main(int argc, char **argv) {
 			}
 			break;
 
+		case 'k':
+			if (strcmp(optarg, "g") == 0) {
+				cluster_maxzoom = MAX_ZOOM - 1;
+				guess_cluster_maxzoom = true;
+			} else {
+				cluster_maxzoom = atoi_require(optarg, "Cluster maxzoom");
+			}
+			break;
+
 		case 'd':
 			full_detail = atoi_require(optarg, "Full detail");
 			if (full_detail > 30) {
@@ -3180,7 +3197,7 @@ int main(int argc, char **argv) {
 
 	long long file_bbox[4] = {UINT_MAX, UINT_MAX, 0, 0};
 
-	ret = read_input(sources, name ? name : out_mbtiles ? out_mbtiles : out_dir, maxzoom, minzoom, basezoom, basezoom_marker_width, outdb, out_dir, &exclude, &include, exclude_all, filter, droprate, buffer, tmpdir, gamma, read_parallel, forcetable, attribution, gamma != 0, file_bbox, prefilter, postfilter, description, guess_maxzoom, &attribute_types, argv[0], &attribute_accum, attribute_descriptions, commandline);
+	ret = read_input(sources, name ? name : out_mbtiles ? out_mbtiles : out_dir, maxzoom, minzoom, basezoom, basezoom_marker_width, outdb, out_dir, &exclude, &include, exclude_all, filter, droprate, buffer, tmpdir, gamma, read_parallel, forcetable, attribution, gamma != 0, file_bbox, prefilter, postfilter, description, guess_maxzoom, guess_cluster_maxzoom, &attribute_types, argv[0], &attribute_accum, attribute_descriptions, commandline);
 
 	if (outdb != NULL) {
 		mbtiles_close(outdb, argv[0]);
diff --git a/main.hpp b/main.hpp
index 8d89ab63c..058c43418 100644
--- a/main.hpp
+++ b/main.hpp
@@ -47,6 +47,7 @@ extern size_t TEMP_FILES;
 extern size_t max_tile_size;
 extern size_t max_tile_features;
 extern int cluster_distance;
+extern int cluster_maxzoom;
 extern std::string attribute_for_id;
 
 int mkstemp_cloexec(char *name);
diff --git a/tile.cpp b/tile.cpp
index cf419825a..6f70fa2ba 100644
--- a/tile.cpp
+++ b/tile.cpp
@@ -1881,7 +1881,7 @@ long long write_tile(FILE *geoms, std::atomic<long long> *geompos_in, char *meta
 				}
 			}
 
-			if (additional[A_CLUSTER_DENSEST_AS_NEEDED] || cluster_distance != 0) {
+			if (z <= cluster_maxzoom && (additional[A_CLUSTER_DENSEST_AS_NEEDED] || cluster_distance != 0)) {
 				indices.push_back(sf.index);
 				if ((sf.index < merge_previndex || sf.index - merge_previndex < mingap) && find_partial(partials, sf, which_partial, layer_unmaps)) {
 					partials[which_partial].clustered++;