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. For the latest version,
# see: https://gist.github.com/bbest/7650602.
#
# This generates the following shapefiles in the output directory:
#   rgns_[offshore|inland]{distance}{units}_[gcs|mol]
# where:
#   rgns = regions
#   offshore = zone extending from shore to the EEZ
#   inland = zone extending from shore inland
#   distance = integer of distance {optional}
#   units = one of: nm (nautical miles), km (kilometers), mi (miles) {optional}
#   gcs = geographic coordinate system
#   mol = Mollweide projection, used to generate Thiessen polygons within a viable global projection
#
# If distance and units are not specified, then the regions are for the full extent,
# ie the entirety of that state's inland area or offshore area extending out to the EEZ.
#
# All shapefiles contain the region id (rgn_id) and name (rgn_name).
# The region id is automatically generated by ascending latitude coordinate in Mollweide projection.
# The rgns_*_gcs shapefiles also have a geodesic area calculated (area_km2).
# 
# Run on cmd: C:\Python27\ArcGISx6410.2\python.exe N:\model\CN-NCEAS-Regions\model.py

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

# set your own paths
lwd      = 'E:/bbest/CN-NCEAS-Regions'    # local working directory
gdb      = lwd+'/scratch.gdb'             # scratch file geodatabase (LOCAL)
rwd      = 'N:/model/CN-NCEAS-Regions'
eez      = rwd+'/tmp/eez_v7_gcs.shp'      # Exclusive Economic Zones (http://marineregions.org)
eezland  = rwd+'/tmp/EEZ_land_v1.shp'     # EEZ plus land (http://marineregions.org)
gadm     = rwd+'/tmp/gadm2.gdb/gadm2'     # Global Administrative Areas (http://gadm.org)
outdir   = rwd+'/data'                    # output directory

# set these variables
country = 'China'
buffers = ['offshore3nm','offshore100nm','offshore1km','inland1km','inland25km']

# buffer units dictionary
buf_units_d = {'nm':'NauticalMiles',
               'km':'Kilometers',
               'mi':'Miles'}

# 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')

# get regions out to eez as full regions offshore
print 'eez_mol...'
arcpy.Intersect_analysis(['eez_mol', 'thiessen_mol'], 'thiessen_eez_inx', 'NO_FID')
arcpy.Dissolve_management('thiessen_eez_inx', 'rgns_offshore_mol', 'NAME_1')

# rgns_offshore: rename NAME_1 to rgn_name
print 'rgns_offshore'
arcpy.AddField_management('rgns_offshore_mol', 'rgn_name', 'TEXT')
arcpy.CalculateField_management('rgns_offshore_mol', 'rgn_name', '!NAME_1!', 'PYTHON_9.3')
arcpy.DeleteField_management('rgns_offshore_mol', 'NAME_1')

# rgns_offshore: assign rgn_id by ascending y coordinate
arcpy.AddField_management('rgns_offshore_mol', 'centroid_y', 'FLOAT')
arcpy.CalculateField_management('rgns_offshore_mol', 'centroid_y', '!shape.centroid.y!', 'PYTHON_9.3')
a = arcpy.da.TableToNumPyArray('rgns_offshore_mol', ['centroid_y','rgn_name'])
a.sort(order=['centroid_y'], axis=0)
a = recfunctions.append_fields(a, 'rgn_id', range(1, a.size+1), usemask=False)
arcpy.da.ExtendTable('rgns_offshore_mol', 'rgn_name', a[['rgn_name','rgn_id']], 'rgn_name', append_only=False)
arcpy.DeleteField_management('rgns_offshore_mol', 'centroid_y')

# rgns_inland
print 'rgns_inland'
arcpy.MakeFeatureLayer_management('states_mol', 'lyr_states_mol')
arcpy.SelectLayerByLocation_management('lyr_states_mol', 'BOUNDARY_TOUCHES', 'eez_mol')
arcpy.CopyFeatures_management('lyr_states_mol', 'rgns_inland_mol')

# rgns_inland: rename NAME_1 to rgn_name
arcpy.AddField_management('rgns_inland_mol', 'rgn_name', 'TEXT')
arcpy.CalculateField_management('rgns_inland_mol', 'rgn_name', '!NAME_1!', 'PYTHON_9.3')
arcpy.DeleteField_management('rgns_inland_mol', 'NAME_1')

# rgns_inland: assign rgn_id
arcpy.da.ExtendTable('rgns_inland_mol', 'rgn_name', a[['rgn_name','rgn_id']], 'rgn_name', append_only=False)

# buffer
for buf in buffers:
    print buf
    rgns_buf_mol = 'rgns_%s_mol' % buf
    buf_zone, buf_dist, buf_units = re.search('(\\D+)(\\d+)(\\D+)', buf).groups()    
    
    if (buf_zone == 'offshore'):
        arcpy.Buffer_analysis('rgns_inland_mol'  , 'buf', '%s %s' % (buf_dist, buf_units_d[buf_units]), 'FULL', 'ROUND', 'ALL')
        arcpy.Intersect_analysis(['buf', 'rgns_offshore_mol'], 'buf_inx', 'NO_FID')
    elif (buf_zone == 'inland'):
        arcpy.Buffer_analysis('rgns_offshore_mol', 'buf', '%s %s' % (buf_dist, buf_units_d[buf_units]), 'FULL', 'ROUND', 'ALL')
        arcpy.Intersect_analysis(['buf', 'rgns_inland_mol'  ], 'buf_inx', 'NO_FID')
    else:
        stop('The buf_zone "%s" is not handled by this function.' % buf_zone)        
    arcpy.Dissolve_management('buf_inx', rgns_buf_mol, ['rgn_id','rgn_name'])
                                                            
# project and calculate area of rgns_*_mol
arcpy.RefreshCatalog(gdb)
for fc in sorted(arcpy.ListFeatureClasses('rgns_*_mol')):
    arcpy.Project_management(fc, fc.replace('_mol', '_gcs'), sr_gcs)
    arcpy.AddField_management(fc, 'area_km2', 'FLOAT')
    arcpy.CalculateField_management(fc, 'area_km2', '!shape.geodesicArea@squarekilometers!', 'PYTHON_9.3')

# copy to shapefiles of rgns_*
arcpy.RefreshCatalog(gdb)
for fc in sorted(arcpy.ListFeatureClasses('rgns_*')):
    if not arcpy.Exists(fc_shp): arcpy.CopyFeatures_management(fc, '%s/%s.shp' % (outdir, fc))