In ESRI geoprocessing terminology, there is an overlay function called identity. This is a very useful function to call when you want to keep all the original geometry boundaries of ONLY the input features combined with an intersection of input features.
This boils down to a formula that calls for both difference and intersect. We first find the difference (input feature - intersection), then add the intersection to create our results as follows:
(input feature – intersection) + intersection = result
For all you curious folks who want to learn how to do this, type out the following code; it will help your muscle memory:
##!/usr/bin/env python # -*- coding: utf-8 -*- from shapely.geometry import asShape, MultiPolygon from utils import shp2_geojson_obj, out_geoj, write_wkt from os.path import realpath def create_polys(shp_data): """ :param shp_data: input GeoJSON :return: MultiPolygon Shapely...