diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index 136575a95..56c12ebae 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -2093,11 +2093,10 @@ def _write_supply_demand_vector(source_aoi_vector_path, feature_attrs,
     Returns:
         ``None``
     """
-    source_vector = ogr.Open(source_aoi_vector_path)
-    driver = ogr.GetDriverByName('GPKG')
-    driver.CopyDataSource(source_vector, target_aoi_vector_path)
-    source_vector = None
-    driver = None
+    gdal.VectorTranslate(
+        target_aoi_vector_path, source_aoi_vector_path,
+        format='GPKG',
+        preserveFID=True)
 
     target_vector = gdal.OpenEx(target_aoi_vector_path, gdal.GA_Update)
     target_layer = target_vector.GetLayer()
diff --git a/src/natcap/invest/wind_energy.py b/src/natcap/invest/wind_energy.py
index 252e4811d..d9e465157 100644
--- a/src/natcap/invest/wind_energy.py
+++ b/src/natcap/invest/wind_energy.py
@@ -1758,10 +1758,9 @@ def _calculate_land_to_grid_distance(
     # Copy the point vector
     _, driver_name = _get_file_ext_and_driver_name(
         target_land_vector_path)
-    base_land_vector = ogr.Open(base_land_vector_path, gdal.OF_VECTOR)
-    driver = ogr.GetDriverByName(driver_name)
-    driver.CopyDataSource(base_land_vector, target_land_vector_path)
-    base_land_vector = None
+    gdal.VectorTranslate(
+        target_land_vector_path, base_land_vector_path,
+        format=driver_name)
 
     target_land_vector = gdal.OpenEx(
         target_land_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
diff --git a/tests/test_coastal_vulnerability.py b/tests/test_coastal_vulnerability.py
index a6e6f6088..3ca634c75 100644
--- a/tests/test_coastal_vulnerability.py
+++ b/tests/test_coastal_vulnerability.py
@@ -931,14 +931,8 @@ class CoastalVulnerabilityTests(unittest.TestCase):
             REGRESSION_DATA, 'coastal_exposure.gpkg')
 
         # This input gets modified in place, so first copy to working dir
-        # I'm using GPKG driver to copy because that driver may have problems
-        # updating a file created by a different GPKG driver version, and the
-        # version used is dependent on GDAL version.
-        # https://gdal.org/drivers/vector/gpkg.html
-        base_shore_point_vector = ogr.Open(base_vector_path)
-        gpkg_driver = ogr.GetDriverByName('GPKG')
-        gpkg_driver.CopyDataSource(
-            base_shore_point_vector, target_point_vector_path)
+        base_shore_point_vector = gdal.OpenEx(base_vector_path, gdal.OF_VECTOR)
+        gdal.VectorTranslate(target_point_vector_path, base_shore_point_vector)
 
         coastal_vulnerability.calculate_final_risk(
             target_point_vector_path, target_point_csv_path)
@@ -968,17 +962,9 @@ class CoastalVulnerabilityTests(unittest.TestCase):
         # This gpkg has a feature with an empty field value for 'R_slr'
         # The function modifies the file in place, so copy to test workspace
         # first.
-
-        # I'm using GPKG driver to copy because that driver may have problems
-        # updating a file created by a different GPKG driver version, and the
-        # version used is dependent on GDAL version.
-        # https://gdal.org/drivers/vector/gpkg.html
         base_vector_path = os.path.join(
             REGRESSION_DATA, 'test_missing_values.gpkg')
-        base_shore_point_vector = ogr.Open(base_vector_path)
-        gpkg_driver = ogr.GetDriverByName('GPKG')
-        gpkg_driver.CopyDataSource(
-            base_shore_point_vector, target_vector_path)
+        gdal.VectorTranslate(target_vector_path, base_vector_path)
 
         coastal_vulnerability.calculate_final_risk(
             target_vector_path, target_csv_path)
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index a3d3d7b43..881e354c9 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -904,8 +904,14 @@ class UNATests(unittest.TestCase):
         from natcap.invest import urban_nature_access
         args = _build_model_args(self.workspace_dir)
 
+        admin_vector = gdal.OpenEx(args['admin_boundaries_vector_path'])
+        admin_layer = admin_vector.GetLayer()
+        fid = admin_layer.GetNextFeature().GetFID()
+        admin_layer = None
+        admin_vector = None
+
         feature_attrs = {
-            0: {
+            fid: {
                 'my-field-1': float(1.2345),
                 'my-field-2': numpy.float32(2.34567),
                 'my-field-3': numpy.float64(3.45678),
@@ -924,10 +930,10 @@ class UNATests(unittest.TestCase):
             vector = gdal.OpenEx(target_vector_path)
             self.assertEqual(vector.GetLayerCount(), 1)
             layer = vector.GetLayer()
-            self.assertEqual(len(layer.schema), len(feature_attrs[0]))
+            self.assertEqual(len(layer.schema), len(feature_attrs[fid]))
             self.assertEqual(layer.GetFeatureCount(), 1)
-            feature = layer.GetFeature(0)
-            for field_name, expected_field_value in feature_attrs[0].items():
+            feature = layer.GetFeature(fid)
+            for field_name, expected_field_value in feature_attrs[fid].items():
                 self.assertEqual(
                     feature.GetField(field_name), expected_field_value)
         finally: