diff --git a/StreamCat.py b/StreamCat.py
index 80e4469..436b172 100644
--- a/StreamCat.py
+++ b/StreamCat.py
@@ -44,6 +44,7 @@
     OUT_DIR,
     PCT_FULL_FILE,
     PCT_FULL_FILE_RP100,
+    USER_ZONES,
 )
 from StreamCat_functions import (
     Accumulation,
@@ -71,7 +72,7 @@
 
 if not os.path.exists(ACCUM_DIR):
     # TODO: work out children OR bastards only
-    makeNumpyVectors(inter_vpu, NHD_DIR)
+    makeNumpyVectors(inter_vpu, NHD_DIR, USER_ZONES)
 
 INPUTS = np.load(ACCUM_DIR +"/vpu_inputs.npy", allow_pickle=True).item()
 
diff --git a/StreamCat_functions.py b/StreamCat_functions.py
index c714c40..a1e1d8a 100644
--- a/StreamCat_functions.py
+++ b/StreamCat_functions.py
@@ -1273,7 +1273,7 @@ def make_all_cat_comids(nhd, inputs):
     return set(all_comids)  # RETURN A SET!
 
 
-def makeNumpyVectors(inter_tbl, nhd):
+def makeNumpyVectors(inter_tbl, nhd, user_zones):
     """
     Uses the NHD tables to create arrays of upstream catchments which are used
     in the Accumulation function
@@ -1284,7 +1284,7 @@ def makeNumpyVectors(inter_tbl, nhd):
     nhd         : directory where NHD is stored
     """
     os.mkdir("accum_npy")
-    inputs = nhd_dict(nhd)
+    inputs = nhd_dict(nhd, user_zones=user_zones)
     all_comids = make_all_cat_comids(nhd, inputs)
     print("Making numpy files in zone...", end="", flush=True)
     for zone, hr in inputs.items():
@@ -1337,7 +1337,7 @@ def makeNumpyVectors(inter_tbl, nhd):
 ##############################################################################
 
 
-def nhd_dict(nhd, unit="VPU"):
+def nhd_dict(nhd, unit="VPU", user_zones=None):
     """
     __author__ =  "Rick Debbout <debbout.rick@epa.gov>"
     Creates an OrderdDict for looping through regions of the NHD to carry
@@ -1354,7 +1354,12 @@ def nhd_dict(nhd, unit="VPU"):
     """
 
     inputs = OrderedDict()
+    if user_zones:  # Use user specified zones
+        inputs |= user_zones
+        np.save("./accum_npy/vpu_inputs.npy", inputs)
+        return inputs
     bounds = dbf2DF(f"{nhd}/NHDPlusGlobalData/BoundaryUnit.dbf")
+    # Drop Hawaii and Cayman Islands.
     remove = bounds.loc[bounds.DRAINAGEID.isin(["HI", "CI"])].index
     bounds = bounds.drop(remove, axis=0)
     if unit == "VPU":
@@ -1406,7 +1411,9 @@ def findUpstreamNpy(zone, com, numpy_dir):
 
 
 def dbf2DF(f, upper=True):
-    data = gpd.read_file(f).drop("geometry", axis=1)
+    data = gpd.read_file(f)
+    if "geometry" in data:
+        data.drop("geometry", axis=1, inplace=True)
     if upper is True:
         data.columns = data.columns.str.upper()
     return data
diff --git a/stream_cat_config.py.template b/stream_cat_config.py.template
index 3dbab08..cff97b9 100644
--- a/stream_cat_config.py.template
+++ b/stream_cat_config.py.template
@@ -18,6 +18,9 @@ STATES_FILE = "/path/to/file/tl_2008_us_state.shp"
 
 ACCUM_DIR = "path/to/local/repository/accump_npy/"
 
+# to run other than all NHD zones, set this dict to e.g. {"04": "GL", "12": "TX"}
+# keys are UnitID and values are DrainageID, see ...\NHDPlusGlobalData\BoundaryUnit.dbf
+USER_ZONES = {}
 
 # location to write out accumulated StreamCat data <- this is intermediate
 OUT_DIR = "/path/to/write/out/files/to"