Feature request: Filter features by intersection with another layer

[mhkeller]

Let's say you have two layers, one congressional district in a state and every zip code in that state. I want to be able to create a layer of zip codes that overlap or touch the congressional district. This would be similar to PostGIS's ST_Intersects.

It could be its own function such as -intersects other_layer target=layer or could be rolled into -filter such as -filter where=intersects(other_layer) target=layer. If this exists already in some capacity, let me know—I wasn't able to find.

The behavior is different from -clip in that there isn't any change to the geometry. In the zip code example above, if you have one zip code straddling the congressional boundary, filtering by intersection would return the zip code feature in its entirety, not just the portion of it clipped to the boundary, as in -clip.

[mbloch]

Mapshaper doesn't have this command currently, but it sounds useful. I'll put it on my list. Thanks for the suggestion :)

[d5xtgr]

Here's a diff that adds partial support for this. I built it on the existing clip/erase logic, so it only works if the test layer is a polygon layer. This is enough for your example, and my separate use case, but not, for instance, finding the countries through which a river (polyline) flows, nor finding what country contains a certain city (point).

To do this properly and support polyline and point test layers is more than I'm presently capable of with my understanding of the codebase. If no one else takes up that challenge, I might attempt it in future once I'm a bit more comfortable with how things work internally.

--- /usr/lib/node_modules/mapshaper/mapshaper.js	2018-09-12 20:09:46.000000000 -0400
+++ ./mapshaper.js	2018-09-13 16:15:12.761816165 -0400
@@ -9784,7 +9784,7 @@
   // Need to expose clip/erase routes in both directions by setting route
   // in both directions to visible -- this is how cut-out shapes are detected
   // Or-ing with 0x11 makes both directions visible (so reverse paths will block)
-  internal.openArcRoutes(clipShapes, arcs, clipFlags, type == 'clip', type == 'erase', !!"dissolve", 0x11);
+  internal.openArcRoutes(clipShapes, arcs, clipFlags, type == 'clip' || type == 'overlap', type == 'erase', !!"dissolve", 0x11);
 
   var index = new PathIndex(clipShapes, arcs);
   var clippedShapes = targetShapes.map(function(shape, i) {
@@ -9805,7 +9805,7 @@
   targetShapes.forEach(function(shape, shapeId) {
     var paths = shape ? findInteriorPaths(shape, type, index) : null;
     if (paths) {
-      clippedShapes[shapeId] = (clippedShapes[shapeId] || []).concat(paths);
+      clippedShapes[shapeId] = (clippedShapes[shapeId] || []).concat((type == 'overlap') ? shape : paths);
     }
   });
 
@@ -9814,6 +9814,7 @@
   function clipPolygon(shape, type, index) {
     var dividedShape = [],
         clipping = type == 'clip',
+        overlapping = type == 'overlap',
         erasing = type == 'erase';
 
     // open pathways for entire polygon rather than one ring at a time --
@@ -9833,7 +9834,7 @@
           // check if it is contained (assumes clip shapes are dissolved)
           if (clipArcTouches === 0 || clipArcUses === 0) { //
             var contained = index.pathIsEnclosed(path);
-            if (clipping && contained || erasing && !contained) {
+            if (clipping && contained || overlapping && contained || erasing && !contained) {
               dividedShape.push(path);
             }
             // TODO: Consider breaking if polygon is unchanged
@@ -9852,7 +9853,7 @@
       usedClipArcs = [];
     }
 
-    return dividedShape.length === 0 ? null : dividedShape;
+    return dividedShape.length === 0 ? null : (type == 'overlap' ? shape : dividedShape);
   }
 
   function routeIsActive(id) {
@@ -9998,7 +10000,7 @@
   function clipPolyline(shp) {
     var clipped = null;
     if (shp) clipped = shp.reduce(clipPath, []);
-    return clipped && clipped.length > 0 ? clipped : null;
+    return clipped && clipped.length > 0 ? (type == 'overlap' ? shp : clipped) : null;
   }
 
   function clipPath(memo, path) {
@@ -10007,7 +10009,7 @@
     for (var i=0; i<path.length; i++) {
       arcId = path[i];
       enclosed = index.arcIsEnclosed(arcId);
-      if (enclosed && type == 'clip' || !enclosed && type == 'erase') {
+      if (enclosed && type == 'clip' || enclosed && type == 'overlap' || !enclosed && type == 'erase') {
         if (!clippedPath) {
           memo.push(clippedPath = []);
         }
@@ -10034,7 +10036,7 @@
 
     for (var i=0; i<n; i++) {
       enclosed = index.findEnclosingShape(feat[i]) > -1;
-      if (type == 'clip' && enclosed || type == 'erase' && !enclosed) {
+      if (type == 'clip' && enclosed || type == 'overlap' && enclosed || type == 'erase' && !enclosed) {
         feat2.push(feat[i].concat());
       }
     }
@@ -10535,6 +10537,11 @@
   return internal.clipLayers(target, src, dataset, "erase", opts);
 };
 
+api.overlapLayers = function(target, src, dataset, opts) {
+  return internal.clipLayers(target, src, dataset, "overlap", opts);
+};
+
+
 api.clipLayer = function(targetLyr, src, dataset, opts) {
   return api.clipLayers([targetLyr], src, dataset, opts)[0];
 };
@@ -10543,6 +10550,10 @@
   return api.eraseLayers([targetLyr], src, dataset, opts)[0];
 };
 
+api.overlapLayer = function(targetLyr, src, dataset, opts) {
+  return api.overlapLayers([targetLyr], src, dataset, opts)[0];
+};
+
 api.sliceLayers = function(target, src, dataset, opts) {
   return internal.clipLayers(target, src, dataset, "slice", opts);
 };
@@ -21013,6 +21021,9 @@
       }
       return internal.writeFiles(outputFiles, opts, done);
 
+    } else if (name == 'overlap') {
+      outputLayers = api.overlapLayers(targetLayers, source, targetDataset, opts);
+
     } else if (name == 'point-grid') {
       outputLayers = [api.pointGrid(targetDataset, opts)];
       if (!targetDataset) {
@@ -22409,6 +22420,20 @@
     .option("debug", {type: "flag"})
     .option("target", targetOpt);
 
+  parser.command("overlap")
+    .describe("remove features that do not overlap another layer")
+    .example("$ mapshaper postcodes.shp -overlap ridings.shp -o overlapped.shp")
+    .validate(validateClipOpts)
+    .option("source", {
+      DEFAULT: true,
+      describe: "file or layer containing overlap polygons"
+    })
+    .option("bbox", bboxOpt)
+    .option("name", nameOpt)
+    .option("no-replace", noReplaceOpt)
+    .option("no-snap", noSnapOpt)
+    .option("target", targetOpt);
+
   parser.command("point-grid")
     .describe("create a rectangular grid of points")
     .validate(validateGridOpts)