Performing an identity function (difference + intersection)

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

How to do it...

1. 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...