bbest · December 29, 2015 09:29 · Feb 20, 2014 · Jan 10, 2014 · Nov 27, 2013 · Nov 25, 2013
diff --git a/create_regions.py b/create_regions.py
@@ -34,7 +34,6 @@
 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)
-gshhs    = rwd+'/tmp/GSHHS_f_L1.shp'      # Global Self-consistent, Hierarchical, High-resolution Geography Database (http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html)
 outdir   = rwd+'/data'                    # output directory
 
 # set these variables

diff --git a/create_regions.py b/create_regions.py
@@ -1,23 +1,50 @@
 # 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/2314628.
+# 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
+import arcpy, os, re, numpy
+from numpy.lib import recfunctions
 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      = 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)
+# 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)
+gshhs    = rwd+'/tmp/GSHHS_f_L1.shp'      # Global Self-consistent, Hierarchical, High-resolution Geography Database (http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html)
+outdir   = rwd+'/data'                    # output directory
 
-# variables
-country      = 'Malaysia'
-buffers      = '3;12;50;100'
-buffer_units = 'NauticalMiles'
+# 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)
@@ -30,7 +57,7 @@
 arcpy.env.outputCoordinateSystem = sr_mol
 
 # select
-arcpy.Select_analysis(eez,     'eez_mol'    , '"Country" = \'%s\'' % country)
+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)
 
@@ -61,15 +88,64 @@
 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
-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')
+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')
 
-# cleanup: project *_mol to *_gcs, and delete the rest
+# copy to shapefiles of rgns_*
 arcpy.RefreshCatalog(gdb)
-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)
+for fc in sorted(arcpy.ListFeatureClasses('rgns_*')):
+    if not arcpy.Exists(fc_shp): arcpy.CopyFeatures_management(fc, '%s/%s.shp' % (outdir, fc))
diff --git a/create_regions.py b/create_regions.py
@@ -1,5 +1,6 @@
 # create_regions.py: generate non-overlapping regions extending administrative areas into the ocean 
-# within the exclusive economic zone by a given buffer
+# within the exclusive economic zone by a given buffer. For the latest version, 
+# see: https://gist.github.com/2314628.
 
 # modules
 import arcpy, os

diff --git a/create_regions.py b/create_regions.py
@@ -7,7 +7,6 @@
 
 # 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)
@@ -18,7 +17,6 @@
 country      = 'Malaysia'
 buffers      = '3;12;50;100'
 buffer_units = 'NauticalMiles'
-#bbox_buffer = '100 Kilometers'
 
 # projections
 sr_mol = arcpy.SpatialReference('Mollweide (world)') # projected Mollweide (54009)
@@ -68,6 +66,7 @@
 arcpy.Intersect_analysis(['buffers_eez', 'thiessen_mol'], 'buffers_mol', 'NO_FID')
 
 # cleanup: project *_mol to *_gcs, and delete the rest
+arcpy.RefreshCatalog(gdb)
 for fc in sorted(arcpy.ListFeatureClasses()):
     if fc.endswith('_mol'):
         arcpy.Project_management(fc, fc.replace('_mol', '_gcs'), sr_gcs)

diff --git a/create_regions.py b/create_regions.py
@@ -1,4 +1,5 @@
-# create_regions.py: generate non-overlapping regions extending administrative areas into the ocean within the exclusive economic zone by a given buffer
+# 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

diff --git a/create_regions.py b/create_regions.py
@@ -1,3 +1,5 @@
+# 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

diff --git a/create_regions.py b/create_regions.py
@@ -0,0 +1,72 @@
+# 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)