diff --git a/clip.cpp b/clip.cpp
index ca0af8e3..a6d77af8 100644
--- a/clip.cpp
+++ b/clip.cpp
@@ -1683,7 +1683,7 @@ static bool feature_out(std::vector<tile_feature> const &features, mvt_layer &ou
 			}
 		}
 
-		outlayer.features.push_back(std::move(outfeature));
+		outlayer.features.push_back(std::make_shared<mvt_feature>(outfeature));
 		return true;
 	}
 
@@ -1694,6 +1694,9 @@ static struct preservecmp {
 	bool operator()(const mvt_feature &a, const mvt_feature &b) {
 		return a.seq < b.seq;
 	}
+	bool operator()(const std::shared_ptr<mvt_feature> &a, const std::shared_ptr<mvt_feature> &b) {
+		return a->seq < b->seq;
+	}
 } preservecmp;
 
 struct index_event {
@@ -1819,49 +1822,10 @@ mvt_tile assign_to_bins(mvt_tile &features,
 			std::set<std::string> exclude,
 			std::vector<std::string> exclude_prefix,
 			int buffer, std::vector<clipbbox> const &clipbboxes) {
-	std::vector<index_event> events;
-	key_pool key_pool;
-
 	if (bins.size() == 0) {
 		return mvt_tile();
 	}
 
-	// Index bins
-	for (size_t i = 0; i < bins.size(); i++) {
-		for (size_t j = 0; j < bins[i].features.size(); j++) {
-			long long xmin, ymin, xmax, ymax;
-			unsigned long long start, end;
-
-			get_bbox(bins[i].features[j].geometry, &xmin, &ymin, &xmax, &ymax, z, x, y, bins[i].detail());
-			get_quadkey_bounds(xmin, ymin, xmax, ymax, &start, &end);
-			events.emplace_back(start, index_event::ENTER, i, j, xmin, ymin, xmax, ymax);
-			events.emplace_back(end, index_event::EXIT, i, j, xmin, ymin, xmax, ymax);
-		}
-	}
-
-	std::map<unsigned long long, std::pair<size_t, size_t>> fid_to_feature;
-
-	// Index points
-	for (size_t i = 0; i < features.layers.size(); i++) {
-		for (size_t j = 0; j < features.layers[i].features.size(); j++) {
-			long long xmin, ymin, xmax, ymax;
-			unsigned long long start, end;
-
-			if (features.layers[i].features[j].geometry.size() > 0) {
-				if (features.layers[i].features[j].has_id) {
-					fid_to_feature.emplace(features.layers[i].features[j].id, std::make_pair(i, j));
-				}
-
-				get_bbox(features.layers[i].features[j].geometry, &xmin, &ymin, &xmax, &ymax, z, x, y, features.layers[i].detail());
-				get_quadkey_bounds(xmin, ymin, xmax, ymax, &start, &end);
-				events.emplace_back(start, index_event::CHECK, i, j, xmin, ymin, xmax, ymax);
-			}
-		}
-	}
-
-	std::sort(events.begin(), events.end());
-	std::set<active_bin> active;
-
 	mvt_layer outlayer;
 	outlayer.extent = bins[0].extent;
 	outlayer.version = 2;
@@ -1869,40 +1833,43 @@ mvt_tile assign_to_bins(mvt_tile &features,
 
 	std::vector<std::vector<tile_feature>> outfeatures;
 
-	for (auto &e : events) {
-		if (e.kind == index_event::ENTER) {
-			active_bin a(e.layer, e.feature);
-			a.xmin = e.xmin;
-			a.ymin = e.ymin;
-			a.xmax = e.xmax;
-			a.ymax = e.ymax;
+	if (bin_by_id_list.size() > 0) {
+		std::map<unsigned long long, std::pair<size_t, size_t>> fid_to_feature;
+		for (size_t i = 0; i < features.layers.size(); i++) {
+			for (size_t j = 0; j < features.layers[i].features.size(); j++) {
+				if (features.layers[i].features[j]->geometry.size() > 0) {
+					if (features.layers[i].features[j]->has_id) {
+						fid_to_feature.emplace(features.layers[i].features[j]->id, std::make_pair(i, j));
+					}
+				}
+			}
+		}
 
-			const mvt_feature &bin = bins[e.layer].features[e.feature];
+		for (auto &binlayer : bins) {
+			size_t seq = 0;
+			for (auto &bin : binlayer.features) {
+				{
+					tile_feature outfeature;
+					for (auto const &g : bin->geometry) {
+						outfeature.geom.emplace_back(g.op, g.x, g.y);
+					}
+					outfeature.t = bin->type;
+					outfeature.has_id = bin->has_id;
+					outfeature.id = bin->id;
+					outfeature.tags = bin->tags;
+					outfeature.layer = &binlayer;
+					outfeature.seq = seq++;
 
-			{
-				tile_feature outfeature;
-				for (auto const &g : bin.geometry) {
-					outfeature.geom.emplace_back(g.op, g.x, g.y);
+					outfeatures.push_back({std::move(outfeature)});
 				}
-				outfeature.t = bin.type;
-				outfeature.has_id = bin.has_id;
-				outfeature.id = bin.id;
-				outfeature.tags = bin.tags;
-				outfeature.layer = &bins[e.layer];
-				outfeature.seq = e.feature;
 
-				a.outfeature = outfeatures.size();
-				outfeatures.push_back({std::move(outfeature)});
-			}
-
-			if (bin_by_id_list.size() > 0) {
-				for (size_t k = 0; k < bin.tags.size(); k += 2) {
-					if (bins[e.layer].keys[bin.tags[k]] == bin_by_id_list) {
-						std::vector<size_t> ids = parse_ids_string(bins[e.layer].values[bin.tags[k + 1]]);
+				for (size_t k = 0; k < bin->tags.size(); k += 2) {
+					if (binlayer.keys[bin->tags[k]] == bin_by_id_list) {
+						std::vector<size_t> ids = parse_ids_string(binlayer.values[bin->tags[k + 1]]);
 						for (auto &id : ids) {
 							auto f = fid_to_feature.find(id);
 							if (f != fid_to_feature.end()) {
-								mvt_feature &feature = features.layers[f->second.first].features[f->second.second];
+								mvt_feature &feature = *features.layers[f->second.first].features[f->second.second];
 								if (feature.geometry.size() > 0) {
 									tile_feature outfeature;
 									for (auto const &g : feature.geometry) {
@@ -1913,8 +1880,8 @@ mvt_tile assign_to_bins(mvt_tile &features,
 									outfeature.has_id = feature.has_id;
 									outfeature.id = feature.id;
 									outfeature.tags = feature.tags;
-									outfeature.layer = &features.layers[e.layer];
-									outfeature.seq = e.feature;
+									outfeature.layer = &features.layers[f->second.first];
+									outfeature.seq = f->second.second;
 									outfeatures.back().push_back(std::move(outfeature));
 								}
 							}
@@ -1923,67 +1890,128 @@ mvt_tile assign_to_bins(mvt_tile &features,
 					}
 				}
 			}
+		}
+	} else {
+		std::vector<index_event> events;
 
-			active.insert(std::move(a));
-		} else if (e.kind == index_event::CHECK) {
-			if (bin_by_id_list.size() > 0) {
-				continue;  // only bin by id, not geometrically
+		// Index bins
+		for (size_t i = 0; i < bins.size(); i++) {
+			for (size_t j = 0; j < bins[i].features.size(); j++) {
+				long long xmin, ymin, xmax, ymax;
+				unsigned long long start, end;
+
+				get_bbox(bins[i].features[j]->geometry, &xmin, &ymin, &xmax, &ymax, z, x, y, bins[i].detail());
+				get_quadkey_bounds(xmin, ymin, xmax, ymax, &start, &end);
+				events.emplace_back(start, index_event::ENTER, i, j, xmin, ymin, xmax, ymax);
+				events.emplace_back(end, index_event::EXIT, i, j, xmin, ymin, xmax, ymax);
 			}
+		}
 
-			auto const &feature = features.layers[e.layer].features[e.feature];
+		// Index points
+		for (size_t i = 0; i < features.layers.size(); i++) {
+			for (size_t j = 0; j < features.layers[i].features.size(); j++) {
+				long long xmin, ymin, xmax, ymax;
+				unsigned long long start, end;
 
-			if (feature.geometry.size() == 0) {
-				// already assigned by ID
-				continue;
+				if (features.layers[i].features[j]->geometry.size() > 0) {
+					get_bbox(features.layers[i].features[j]->geometry, &xmin, &ymin, &xmax, &ymax, z, x, y, features.layers[i].detail());
+					get_quadkey_bounds(xmin, ymin, xmax, ymax, &start, &end);
+					events.emplace_back(start, index_event::CHECK, i, j, xmin, ymin, xmax, ymax);
+				}
 			}
+		}
 
-			// if we can't find a real match,
-			// assign points to the most nearby bin
-			ssize_t which_outfeature = outfeatures.size() - 1;
+		std::sort(events.begin(), events.end());
+		std::set<active_bin> active;
 
-			for (auto const &a : active) {
-				auto const &bin = bins[a.layer].features[a.feature];
+		for (auto &e : events) {
+			if (e.kind == index_event::ENTER) {
+				active_bin a(e.layer, e.feature);
+				a.xmin = e.xmin;
+				a.ymin = e.ymin;
+				a.xmax = e.xmax;
+				a.ymax = e.ymax;
 
-				if (bbox_intersects(e.xmin, e.ymin, e.xmax, e.ymax,
-						    a.xmin, a.ymin, a.xmax, a.ymax)) {
-					if (pnpoly_mp(bin.geometry, feature.geometry[0].x, feature.geometry[0].y)) {
-						which_outfeature = a.outfeature;
-						break;
+				const mvt_feature &bin = *bins[e.layer].features[e.feature];
+
+				{
+					tile_feature outfeature;
+					for (auto const &g : bin.geometry) {
+						outfeature.geom.emplace_back(g.op, g.x, g.y);
 					}
+					outfeature.t = bin.type;
+					outfeature.has_id = bin.has_id;
+					outfeature.id = bin.id;
+					outfeature.tags = bin.tags;
+					outfeature.layer = &bins[e.layer];
+					outfeature.seq = e.feature;
+
+					a.outfeature = outfeatures.size();
+					outfeatures.push_back({std::move(outfeature)});
 				}
-			}
 
-			if (which_outfeature >= 0) {
-				tile_feature outfeature;
-				for (auto const &g : feature.geometry) {
-					outfeature.geom.emplace_back(g.op, g.x, g.y);
-				}
-				outfeature.t = feature.type;
-				outfeature.has_id = feature.has_id;
-				outfeature.id = feature.id;
-				outfeature.tags = feature.tags;
-				outfeature.layer = &features.layers[e.layer];
-				outfeature.seq = e.feature;
-				outfeatures[which_outfeature].push_back(std::move(outfeature));
-			}
-		} else /* EXIT */ {
-			auto const &found = active.find({e.layer, e.feature});
-			if (found != active.end()) {
-				active.erase(found);
-			} else {
-				fprintf(stderr, "event mismatch: can't happen\n");
-				exit(EXIT_IMPOSSIBLE);
+				active.insert(std::move(a));
+			} else if (e.kind == index_event::CHECK) {
+				if (bin_by_id_list.size() > 0) {
+					continue;  // only bin by id, not geometrically
+				}
+
+				auto const &feature = *features.layers[e.layer].features[e.feature];
+
+				if (feature.geometry.size() == 0) {
+					// already assigned by ID
+					continue;
+				}
+
+				// if we can't find a real match,
+				// assign points to the most nearby bin
+				ssize_t which_outfeature = outfeatures.size() - 1;
+
+				for (auto const &a : active) {
+					auto const &bin = *bins[a.layer].features[a.feature];
+
+					if (bbox_intersects(e.xmin, e.ymin, e.xmax, e.ymax,
+							    a.xmin, a.ymin, a.xmax, a.ymax)) {
+						if (pnpoly_mp(bin.geometry, feature.geometry[0].x, feature.geometry[0].y)) {
+							which_outfeature = a.outfeature;
+							break;
+						}
+					}
+				}
+
+				if (which_outfeature >= 0) {
+					tile_feature outfeature;
+					for (auto const &g : feature.geometry) {
+						outfeature.geom.emplace_back(g.op, g.x, g.y);
+					}
+					outfeature.t = feature.type;
+					outfeature.has_id = feature.has_id;
+					outfeature.id = feature.id;
+					outfeature.tags = feature.tags;
+					outfeature.layer = &features.layers[e.layer];
+					outfeature.seq = e.feature;
+					outfeatures[which_outfeature].push_back(std::move(outfeature));
+				}
+			} else /* EXIT */ {
+				auto const &found = active.find({e.layer, e.feature});
+				if (found != active.end()) {
+					active.erase(found);
+				} else {
+					fprintf(stderr, "event mismatch: can't happen\n");
+					exit(EXIT_IMPOSSIBLE);
+				}
 			}
 		}
 	}
 
+	key_pool key_pool;
 	for (size_t i = 0; i < outfeatures.size(); i++) {
 		if (outfeatures[i].size() > 1) {
 			if (feature_out(outfeatures[i], outlayer,
 					keep, exclude, exclude_prefix, attribute_accum,
 					accumulate_numeric, key_pool, buffer, true,
 					clipbboxes, z, x, y, NULL)) {
-				mvt_feature &nfeature = outlayer.features.back();
+				mvt_feature &nfeature = *outlayer.features.back();
 				mvt_value val;
 				val.type = mvt_uint;
 				val.numeric_value.uint_value = outfeatures[i].size() - 1;
@@ -2069,7 +2097,8 @@ std::string overzoom(std::vector<source_tile> const &tiles, int nz, int nx, int
 			static const std::string retain_points_multiplier_first = "tippecanoe:retain_points_multiplier_first";
 			static const std::string retain_points_multiplier_sequence = "tippecanoe:retain_points_multiplier_sequence";
 
-			for (auto feature : layer.features) {
+			for (auto feature_ref : layer.features) {
+				auto feature = *feature_ref;  // this needs to be a copy, because tile-join still needs the original unmutated
 				drawvec geom;
 				int t = feature.type;
 
diff --git a/decode.cpp b/decode.cpp
index 534e6219..febdbba9 100644
--- a/decode.cpp
+++ b/decode.cpp
@@ -60,11 +60,11 @@ void do_stats(mvt_tile &tile, size_t size, bool compressed, int z, unsigned x, u
 
 		size_t points = 0, lines = 0, polygons = 0;
 		for (size_t j = 0; j < tile.layers[i].features.size(); j++) {
-			if (tile.layers[i].features[j].type == mvt_point) {
+			if (tile.layers[i].features[j]->type == mvt_point) {
 				points++;
-			} else if (tile.layers[i].features[j].type == mvt_linestring) {
+			} else if (tile.layers[i].features[j]->type == mvt_linestring) {
 				lines++;
-			} else if (tile.layers[i].features[j].type == mvt_polygon) {
+			} else if (tile.layers[i].features[j]->type == mvt_polygon) {
 				polygons++;
 			}
 		}
diff --git a/mvt.cpp b/mvt.cpp
index 16f9f583..3d31df7d 100644
--- a/mvt.cpp
+++ b/mvt.cpp
@@ -278,7 +278,7 @@ bool mvt_tile::decode(const std::string &message, bool &was_compressed) {
 						}
 					}
 
-					layer.features.push_back(std::move(feature));
+					layer.features.push_back(std::make_shared<mvt_feature>(feature));
 					break;
 				}
 
@@ -407,16 +407,16 @@ std::string mvt_tile::encode() {
 			std::string feature_string;
 			protozero::pbf_writer feature_writer(feature_string);
 
-			feature_writer.add_enum(3, layers[i].features[f].type);
+			feature_writer.add_enum(3, layers[i].features[f]->type);
 
-			std::vector<unsigned> sorted_tags = layers[i].features[f].tags;
+			std::vector<unsigned> sorted_tags = layers[i].features[f]->tags;
 			for (size_t v = 1; v < sorted_tags.size(); v += 2) {
 				sorted_tags[v] = mapping[sorted_tags[v]];
 			}
 			feature_writer.add_packed_uint32(2, std::begin(sorted_tags), std::end(sorted_tags));
 
-			if (layers[i].features[f].has_id) {
-				feature_writer.add_uint64(1, layers[i].features[f].id);
+			if (layers[i].features[f]->has_id) {
+				feature_writer.add_uint64(1, layers[i].features[f]->id);
 			}
 
 			std::vector<uint32_t> geometry;
@@ -426,7 +426,7 @@ std::string mvt_tile::encode() {
 			int cmd = -1;
 			int length = 0;
 
-			std::vector<mvt_geometry> &geom = layers[i].features[f].geometry;
+			std::vector<mvt_geometry> &geom = layers[i].features[f]->geometry;
 
 			for (size_t g = 0; g < geom.size(); g++) {
 				int op = geom[g].op;
diff --git a/mvt.hpp b/mvt.hpp
index 8f24fc68..7e52363c 100644
--- a/mvt.hpp
+++ b/mvt.hpp
@@ -204,7 +204,7 @@ struct std::hash<mvt_value> {
 struct mvt_layer {
 	int version = 0;
 	std::string name = "";
-	std::vector<mvt_feature> features{};
+	std::vector<std::shared_ptr<mvt_feature>> features{};
 	std::vector<std::string> keys{};
 	std::vector<mvt_value> values{};
 	long long extent = 0;
diff --git a/plugin.cpp b/plugin.cpp
index cfe06593..b6033111 100644
--- a/plugin.cpp
+++ b/plugin.cpp
@@ -119,17 +119,17 @@ std::vector<mvt_layer> parse_layers(int fd, int z, unsigned x, unsigned y, std::
 		}
 
 		for (auto const &feature : layer.features) {
-			if (feature.type == mvt_point) {
+			if (feature->type == mvt_point) {
 				ts->second.points++;
-			} else if (feature.type == mvt_linestring) {
+			} else if (feature->type == mvt_linestring) {
 				ts->second.lines++;
-			} else if (feature.type == mvt_polygon) {
+			} else if (feature->type == mvt_polygon) {
 				ts->second.polygons++;
 			}
 
-			for (size_t i = 0; i + 1 < feature.tags.size(); i += 2) {
-				const std::string &key = layer.keys[feature.tags[i]];
-				const mvt_value &val = layer.values[feature.tags[i + 1]];
+			for (size_t i = 0; i + 1 < feature->tags.size(); i += 2) {
+				const std::string &key = layer.keys[feature->tags[i]];
+				const mvt_value &val = layer.values[feature->tags[i + 1]];
 
 				// Nulls can be excluded here because this is the postfilter
 				// and it is nearly time to create the vector representation
diff --git a/read_json.cpp b/read_json.cpp
index 63329acb..88ec04bc 100644
--- a/read_json.cpp
+++ b/read_json.cpp
@@ -398,7 +398,7 @@ std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int
 				}
 			}
 
-			l->second.features.push_back(feature);
+			l->second.features.push_back(std::make_shared<mvt_feature>(feature));
 		}
 
 		json_free(j);
diff --git a/tile-join.cpp b/tile-join.cpp
index 2411fe06..052c91c0 100644
--- a/tile-join.cpp
+++ b/tile-join.cpp
@@ -243,9 +243,9 @@ void append_tile(std::string message, int z, unsigned x, unsigned y, std::map<st
 				// scaling up is safe while scaling down requires geometry cleaning.
 
 				for (size_t i = 0; i < outlayer.features.size(); i++) {
-					for (size_t j = 0; j < outlayer.features[i].geometry.size(); j++) {
-						outlayer.features[i].geometry[j].x = outlayer.features[i].geometry[j].x * layer.extent / outlayer.extent;
-						outlayer.features[i].geometry[j].y = outlayer.features[i].geometry[j].y * layer.extent / outlayer.extent;
+					for (size_t j = 0; j < outlayer.features[i]->geometry.size(); j++) {
+						outlayer.features[i]->geometry[j].x = outlayer.features[i]->geometry[j].x * layer.extent / outlayer.extent;
+						outlayer.features[i]->geometry[j].y = outlayer.features[i]->geometry[j].y * layer.extent / outlayer.extent;
 					}
 				}
 
@@ -261,7 +261,7 @@ void append_tile(std::string message, int z, unsigned x, unsigned y, std::map<st
 			join_keys.resize(layer.features.size());
 
 			for (size_t f = 0; f < layer.features.size(); f++) {
-				mvt_feature &feat = layer.features[f];
+				mvt_feature &feat = *layer.features[f];
 				join_keys[f].type = mvt_no_such_key;
 
 				for (size_t t = 0; t + 1 < feat.tags.size(); t += 2) {
@@ -289,7 +289,7 @@ void append_tile(std::string message, int z, unsigned x, unsigned y, std::map<st
 		bool features_added_to_layer = false;
 
 		for (size_t f = 0; f < layer.features.size(); f++) {
-			mvt_feature &feat = layer.features[f];
+			mvt_feature &feat = *layer.features[f];
 
 			std::set<std::string> exclude_attributes;
 			if (filter != NULL && !evaluate(feat, layer, filter, exclude_attributes, z, unidecode_data)) {
@@ -520,7 +520,7 @@ void append_tile(std::string message, int z, unsigned x, unsigned y, std::map<st
 
 				features_added++;
 				features_added_to_layer = true;
-				outlayer.features.push_back(outfeature);
+				outlayer.features.push_back(std::make_shared<mvt_feature>(outfeature));
 
 				if (z < tilestats->second.minzoom) {
 					tilestats->second.minzoom = z;
diff --git a/tile.cpp b/tile.cpp
index 312c289a..f6c4e7b1 100644
--- a/tile.cpp
+++ b/tile.cpp
@@ -1354,7 +1354,7 @@ void *run_prefilter(void *v) {
 		}
 
 		decode_meta(sf, tmp_layer, tmp_feature);
-		tmp_layer.features.push_back(tmp_feature);
+		tmp_layer.features.push_back(std::make_shared<mvt_feature>(tmp_feature));
 
 		layer_to_geojson(tmp_layer, 0, 0, 0, false, true, false, true, sf.index, sf.seq, sf.extent, true, state, 0, std::set<std::string>());
 	}
@@ -2587,7 +2587,7 @@ long long write_tile(decompressor *geoms, std::atomic<long long> *geompos_in, ch
 					add_tilestats(layer.name, z, layermaps, tiling_seg, layer_unmaps, "tippecanoe_feature_density", sv);
 				}
 
-				layer.features.push_back(std::move(feature));
+				layer.features.push_back(std::make_shared<mvt_feature>(feature));
 				layer_features[x] = std::make_shared<serial_feature>();
 			}
 
diff --git a/write_json.cpp b/write_json.cpp
index eae5a2f0..fc11dc56 100644
--- a/write_json.cpp
+++ b/write_json.cpp
@@ -270,7 +270,7 @@ void write_coords(json_writer &state, lonlat const &ll, double scale) {
 
 void layer_to_geojson(mvt_layer const &layer, unsigned z, unsigned x, unsigned y, bool comma, bool name, bool zoom, bool write_dropped, unsigned long long index, long long sequence, long long extent, bool complain, json_writer &state, double scale, std::set<std::string> const &include_attr) {
 	for (size_t f = 0; f < layer.features.size(); f++) {
-		mvt_feature const &feat = layer.features[f];
+		mvt_feature const &feat = *layer.features[f];
 
 		state.json_write_hash();
 		state.json_write_string("type");