generate non-overlapping regions extending administrative areas into the ocean within the exclusive economic zone by a given buffer

create_regions.py

# create_regions.py: generate non-overlapping regions extending administrative areas into the ocean 
# within the exclusive economic zone by a given buffer

# modules
import arcpy, os
arcpy.SetLogHistory(True) # C:\Users\bbest\AppData\Roaming\ESRI\Desktop10.2\ArcToolbox\History

# paths
lwd      = r'E:\bbest\neptune_local_edit_src_toolbox_exercises_create_regions'   # local working directory
#gdb      = r'L:\src\toolbox\exercises\create_regions\scratch.gdb'               # scratch file geodatabase (NETWORK)
gdb      = lwd+r'\scratch.gdb'                                                   # scratch file geodatabase (LOCAL)
eez      = r'N:\stable\GL-VLIZ-EEZs_v7\data\eez_v7_gcs.shp'                      # Exclusive Economic Zones (http://marineregions.org)
eezland  = r'N:\stable\GL-VLIZ-EEZs_v7\data\EEZ_land_v1.shp'                     # EEZ plus land (http://marineregions.org)
gadm     = r'N:\stable\GL-GADM-AdminAreas_v2\data\gadm2.gdb\gadm2'               # Global Administrative Areas (http://gadm.org)
gshhs    = r'N:\stable\GL-NGDC-GSHHG\gshhg-shp-2.2.2\GSHHS_shp\f\GSHHS_f_L1.shp' # Global Self-consistent, Hierarchical, High-resolution Geography Database (http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html)

# variables
country      = 'Malaysia'
buffers      = '3;12;50;100'
buffer_units = 'NauticalMiles'
#bbox_buffer = '100 Kilometers'

# projections
sr_mol = arcpy.SpatialReference('Mollweide (world)') # projected Mollweide (54009)
sr_gcs = arcpy.SpatialReference('WGS 1984')          # geographic coordinate system WGS84 (4326)

# environment
if not arcpy.Exists(gdb): arcpy.CreateFileGDB_management(os.path.dirname(gdb), os.path.basename(gdb))
arcpy.env.workspace       = gdb
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = sr_mol

# select
arcpy.Select_analysis(eez,     'eez_mol'    , '"Country" = \'%s\'' % country)
arcpy.Select_analysis(eezland, 'eezland', '"Country" = \'%s\'' % country)
arcpy.Select_analysis(gadm,    'gadm'   , '"NAME_0"  = \'%s\'' % country)

# get administrative land
arcpy.Erase_analysis('eezland', 'eez_mol', 'land_mol')
arcpy.Clip_analysis('gadm', 'land_mol', 'gadm_land')
arcpy.Dissolve_management('gadm_land', 'states', 'NAME_1')

# create theissen polygons used to split slivers
arcpy.Densify_edit("states", 'DISTANCE', '1 Kilometers')
arcpy.FeatureVerticesToPoints_management('states', 'states_pts', 'ALL')

# delete interior points
arcpy.Dissolve_management('states', 'states_dissolved')
arcpy.MakeFeatureLayer_management('states_pts', 'lyr_states_pts')
arcpy.SelectLayerByLocation_management('lyr_states_pts', 'WITHIN_CLEMENTINI', 'states_dissolved')
arcpy.DeleteFeatures_management('lyr_states_pts')

# generate thiessen polygons of gadm for intersecting with land slivers
arcpy.env.extent = 'eezland'
arcpy.CreateThiessenPolygons_analysis('states_pts', 'states_thiessen', 'ALL')
arcpy.Dissolve_management('states_thiessen', 'thiessen_mol', 'NAME_1')
arcpy.RepairGeometry_management('thiessen_mol')
    
# get slivers, which are land but not identified by gadm, intersect with thiessen so break at junctions
arcpy.Erase_analysis('land_mol', 'gadm_land', 'landnotgadm')
arcpy.Intersect_analysis(["landnotgadm", 'thiessen_mol'], 'slivers')
arcpy.Merge_management(['states', 'slivers'], 'states_slivers')
arcpy.Dissolve_management('states_slivers', 'states_mol', 'NAME_1')

# buffer
arcpy.MultipleRingBuffer_analysis('states_mol', 'buffers', buffers, buffer_units, 'distance', 'ALL', 'OUTSIDE_ONLY')
arcpy.Clip_analysis('buffers', 'eez_mol', 'buffers_eez')
arcpy.Intersect_analysis(['buffers_eez', 'thiessen_mol'], 'buffers_mol', 'NO_FID')

# cleanup: project *_mol to *_gcs, and delete the rest
for fc in sorted(arcpy.ListFeatureClasses()):
    if fc.endswith('_mol'):
        arcpy.Project_management(fc, fc.replace('_mol', '_gcs'), sr_gcs)
    else:
        arcpy.Delete_management(fc)