cleanup; remove vector points beyond raster extents; rework pixel-picking test (#1698)
This commit is contained in:
Megan Nissel
parent
eb511d13a8
commit
5d6affde14
|
|
@ -870,10 +870,7 @@ def execute(args):
|
|||
# since min_distance and max_distance will both be ''
|
||||
except (KeyError, ValueError):
|
||||
LOGGER.info('Distance information not provided')
|
||||
mask_by_distance = False
|
||||
|
||||
else:
|
||||
mask_by_distance = True
|
||||
# Clip and project the land polygon shapefile to AOI
|
||||
LOGGER.info('Clip and project land polygon to AOI')
|
||||
land_poly_proj_vector_path = os.path.join(
|
||||
|
|
@ -916,7 +913,11 @@ def execute(args):
|
|||
|
||||
else:
|
||||
LOGGER.info("AOI argument was not selected")
|
||||
mask_by_distance = False
|
||||
if run_valuation:
|
||||
# Guard against trying to run the Valuation model without an AOI
|
||||
LOGGER.warning("Run Valuation was selected but no AOI was provided. "
|
||||
"Please provide an AOI in order to run Valuation.")
|
||||
run_valuation = False
|
||||
|
||||
# Wind point vector that will be the template for the final
|
||||
# shapefile; does not need to be clipped to AOI
|
||||
|
|
@ -966,8 +967,6 @@ def execute(args):
|
|||
mask_keys.append(_DEPTH_FIELD_NAME)
|
||||
write_vector_dependent_task_list.append(create_depth_mask_task)
|
||||
|
||||
# If running valuation, raster alignment will happen after rasterization
|
||||
# of the Harvested values (in the valuation half of the model)
|
||||
if not run_valuation:
|
||||
# Write Depth [and Distance] mask values to Wind Points Shapefile
|
||||
LOGGER.info("Adding mask values to shapefile")
|
||||
|
|
@ -978,8 +977,7 @@ def execute(args):
|
|||
args=(intermediate_wind_point_vector_path,
|
||||
raster_field_to_vector_list,
|
||||
final_wind_point_vector_path),
|
||||
kwargs={'overwrite': True,
|
||||
'mask_keys': mask_keys,
|
||||
kwargs={'mask_keys': mask_keys,
|
||||
'mask_field': _MASK_FIELD_NAME},
|
||||
target_path_list=[final_wind_point_vector_path],
|
||||
task_name='add_masked_vals_to_wind_vector',
|
||||
|
|
@ -1036,7 +1034,7 @@ def execute(args):
|
|||
task_name='mask_harvested_raster',
|
||||
target_path_list=[harvested_masked_path],
|
||||
dependent_task_list=[rasterize_harvested_task,
|
||||
create_depth_mask_task, create_dist_mask_task])
|
||||
create_dist_mask_task])
|
||||
|
||||
# path for final distance transform used in valuation calculations
|
||||
final_dist_raster_path = os.path.join(
|
||||
|
|
@ -1268,8 +1266,7 @@ def execute(args):
|
|||
args=(intermediate_wind_point_vector_path,
|
||||
raster_field_to_vector_list,
|
||||
final_wind_point_vector_path),
|
||||
kwargs={'overwrite': True,
|
||||
'mask_keys': mask_keys,
|
||||
kwargs={'mask_keys': mask_keys,
|
||||
'mask_field': _MASK_FIELD_NAME},
|
||||
target_path_list=[final_wind_point_vector_path],
|
||||
task_name='add_harv_valuation_to_wind_vector',
|
||||
|
|
@ -1282,44 +1279,43 @@ def execute(args):
|
|||
|
||||
def _index_raster_values_to_point_vector(
|
||||
base_point_vector_path, raster_fieldname_list,
|
||||
target_point_vector_path, overwrite=False,
|
||||
mask_keys=[], mask_field=None):
|
||||
target_point_vector_path, mask_keys=[], mask_field=None):
|
||||
"""Add raster values to vector point feature fields.
|
||||
|
||||
This function does two things:
|
||||
- Creates a copy of the base point vector and updates the copy to include
|
||||
the fields in raster_fieldname_list, the values of which are pixel-picked
|
||||
from their corresponding raster. If `mask_keys` are provided, points that
|
||||
would be NODATA based on the raster(s) corresponding to fields in this list
|
||||
will be deleted from the copy.
|
||||
the fields in ``raster_fieldname_list``, the values of which are
|
||||
pixel-picked from their corresponding raster. If ``mask_keys`` are
|
||||
provided, points that would be NODATA based on the raster(s)
|
||||
corresponding to fields in this list will be deleted from the copy.
|
||||
- Optionally modifies the base vector to include an additional field,
|
||||
`mask_field`, that indicates whether or not a point was removed from the
|
||||
copy based on NODATA values in the rasters corresponding to the fieldnames
|
||||
in `mask_keys`.
|
||||
``mask_field``, that indicates whether or not a point was removed from
|
||||
the copy based on NODATA values in the rasters corresponding to the
|
||||
fieldnames in ``mask_keys``.
|
||||
|
||||
If the base vector contains points that fall beyond the extent of the rasters,
|
||||
those points will be removed from the output vector. If ``mask_field`` is
|
||||
provided, the value will remain NULL for those points.
|
||||
|
||||
Args:
|
||||
base_point_vector_path (str): a path to an OGR point vector file.
|
||||
raster_fieldname_list (list): a list of (raster_path, field_name)
|
||||
tuples. The values of rasters in this list will be added to the
|
||||
vector's features under the associated field name. All rasters must
|
||||
already be aligned. An exception will be raised if a field already
|
||||
exists in the base point vector and overwrite is False.
|
||||
already be aligned.
|
||||
target_point_vector_path (str): a path to a shapefile that has the
|
||||
target field name in addition to the existing fields in the base
|
||||
point vector.
|
||||
overwrite (boolean): defaults to False. If True, if a fieldname passed
|
||||
in raster_fieldname_list already exists in the source shapefile,
|
||||
its value will be overwritten in the copy.
|
||||
mask_keys (list): (optional) a list of string field names that appear in
|
||||
raster_fieldname_list, indicating that the associated raster
|
||||
``raster_fieldname_list``, indicating that the associated raster
|
||||
will be used to mask out vector features. Wherever a value in a mask
|
||||
raster equals NODATA, the vector feature will be deleted from the
|
||||
target point vector.
|
||||
mask_field (str): (optional) a field name to add to the base point vector,
|
||||
where the value indicates whether or not the point was masked out in
|
||||
the target point vector, based on the values of rasters associated
|
||||
with the fields in mask_keys. If `mask_field` is provided but
|
||||
`mask_keys` is not, it will be ignored.
|
||||
with the fields in ``mask_keys``. If ``mask_field`` is provided but
|
||||
``mask_keys`` is not, it will be ignored.
|
||||
|
||||
Returns:
|
||||
None
|
||||
|
|
@ -1364,8 +1360,8 @@ def _index_raster_values_to_point_vector(
|
|||
base_vector = gdal.OpenEx(
|
||||
base_point_vector_path, gdal.OF_VECTOR | gdal.GA_Update)
|
||||
base_layer = base_vector.GetLayer()
|
||||
# Since we're mutating the base vector, throw an error if the
|
||||
# supplied fieldname already exists (don't overwrite)
|
||||
# Since we're mutating the base vector, if the
|
||||
# supplied fieldname already exists, don't overwrite it
|
||||
if base_layer.FindFieldIndex(mask_field, True) != -1:
|
||||
raise ValueError(
|
||||
f"A field called '{mask_field}' already exists in the base "
|
||||
|
|
@ -1389,14 +1385,8 @@ def _index_raster_values_to_point_vector(
|
|||
field_defn.SetPrecision(11)
|
||||
|
||||
field_index = target_layer.FindFieldIndex(field_name, True)
|
||||
# If the field exists and overwrite = False, raise an exception
|
||||
if field_index != -1 and not overwrite:
|
||||
raise ValueError(
|
||||
f"'{field_name}' field already exists in the input shapefile and "
|
||||
"overwrite is set to False. Please rename it or remove it "
|
||||
"from the attribute table.")
|
||||
# If the field doesn't already exist, create it
|
||||
elif field_index == -1:
|
||||
if field_index == -1:
|
||||
target_layer.CreateField(field_defn)
|
||||
|
||||
# Create coordinate transformation from vector to raster, to make sure the
|
||||
|
|
@ -1410,6 +1400,10 @@ def _index_raster_values_to_point_vector(
|
|||
vector_coord_trans = utils.create_coordinate_transformer(
|
||||
vector_sr, raster_sr)
|
||||
|
||||
# We'll check encountered features against this list at the end,
|
||||
# for the purposes of removing any points that weren't encountered
|
||||
all_fids = [feat.GetFID() for feat in target_layer]
|
||||
|
||||
# Initialize an R-Tree indexing object with point geom from base_vector
|
||||
def generator_function():
|
||||
for feat in target_layer:
|
||||
|
|
@ -1423,6 +1417,9 @@ def _index_raster_values_to_point_vector(
|
|||
|
||||
vector_idx = index.Index(generator_function(), interleaved=False)
|
||||
|
||||
# For all the features (points) add the proper raster value
|
||||
encountered_fids = set()
|
||||
|
||||
iterables = [pygeoprocessing.iterblocks((base_raster_path_list[index], 1))
|
||||
for index in range(len(base_raster_path_list))]
|
||||
for block_data in zip(*iterables):
|
||||
|
|
@ -1431,9 +1428,6 @@ def _index_raster_values_to_point_vector(
|
|||
block_matrices = [block_data[index][1]
|
||||
for index in range(len(base_raster_path_list))]
|
||||
|
||||
# For all the features (points) add the proper raster value
|
||||
encountered_fids = set()
|
||||
|
||||
block_min_x = raster_min_x + block_info['xoff'] * pixel_size_x
|
||||
block_max_x = raster_min_x + (
|
||||
block_info['win_xsize'] + block_info['xoff']) * pixel_size_x
|
||||
|
|
@ -1501,6 +1495,12 @@ def _index_raster_values_to_point_vector(
|
|||
target_layer.SetFeature(target_vector_feat)
|
||||
feat = None
|
||||
|
||||
# Finally, if any points fall outside of the raster extents, remove them
|
||||
out_of_bounds_fids = set(all_fids).difference(encountered_fids)
|
||||
if out_of_bounds_fids:
|
||||
for vector_fid in out_of_bounds_fids:
|
||||
target_layer.DeleteFeature(vector_fid)
|
||||
|
||||
target_vector.ExecuteSQL('REPACK ' + target_layer.GetName())
|
||||
target_layer = None
|
||||
target_vector = None
|
||||
|
|
|
|||
|
|
@ -404,100 +404,69 @@ class WindEnergyUnitTests(unittest.TestCase):
|
|||
"""WindEnergy: testing 'index_raster_values_to_point_vector' function."""
|
||||
from natcap.invest import wind_energy
|
||||
|
||||
base_vector_path = os.path.join(self.workspace_dir, 'base_vector.shp')
|
||||
points = [
|
||||
{'long': -69.49, 'lati': 43.89, 'pt_id': 0}, # valid
|
||||
{'long': -69.48, 'lati': 43.89, 'pt_id': 1}, # masked out
|
||||
{'long': -69.49, 'lati': 43.88, 'pt_id': 2}, # masked out
|
||||
{'long': -69.48, 'lati': 43.88, 'pt_id': 3}, # valid
|
||||
]
|
||||
driver = gdal.GetDriverByName('ESRI Shapefile')
|
||||
target_vector = driver.Create(base_vector_path, 0, 0, 0, gdal.GDT_Unknown)
|
||||
|
||||
source_sr = osr.SpatialReference()
|
||||
source_sr.SetWellKnownGeogCS("WGS84") # lati / long
|
||||
output_layer = target_vector.CreateLayer('wind_points', source_sr, ogr.wkbPoint)
|
||||
|
||||
for field in points[0].keys():
|
||||
field_defn = ogr.FieldDefn(field, ogr.OFTReal)
|
||||
output_layer.CreateField(field_defn)
|
||||
|
||||
for point in points:
|
||||
geom = ogr.Geometry(ogr.wkbPoint)
|
||||
geom.AddPoint_2D(point['long'], point['lati'])
|
||||
|
||||
output_feature = ogr.Feature(output_layer.GetLayerDefn())
|
||||
|
||||
for field_name in point.keys():
|
||||
field_index = output_feature.GetFieldIndex(field_name)
|
||||
output_feature.SetField(field_index, point[field_name])
|
||||
|
||||
output_feature.SetGeometryDirectly(geom)
|
||||
output_layer.CreateFeature(output_feature)
|
||||
output_feature = None
|
||||
|
||||
output_layer.SyncToDisk()
|
||||
output_layer = None
|
||||
target_vector = None
|
||||
|
||||
srs = osr.SpatialReference()
|
||||
srs.ImportFromEPSG(32619)
|
||||
srs.ImportFromEPSG(3157)
|
||||
projection_wkt = srs.ExportToWkt()
|
||||
origin = (443723.127327877911739, 4956546.905980412848294)
|
||||
pos_x = origin[0]
|
||||
pos_y = origin[1]
|
||||
|
||||
# When rasterizing values, we'll always be working in projected space
|
||||
base_proj_vector_path = os.path.join(
|
||||
self.workspace_dir, 'base_proj_vector.shp')
|
||||
pygeoprocessing.reproject_vector(
|
||||
base_vector_path, projection_wkt, base_proj_vector_path)
|
||||
# Setup parameters for creating point shapefile
|
||||
fields = {'id': ogr.OFTReal}
|
||||
attrs = [
|
||||
{'id': 0}, {'id': 1}, {'id': 2},
|
||||
{'id': 3}, {'id': 4}, {'id': 5}]
|
||||
valid_ids = [1, 5]
|
||||
|
||||
initial_raster_path = os.path.join(
|
||||
self.workspace_dir, 'empty_raster.tif')
|
||||
# Harvested values are rasterized from the vector;
|
||||
# create our base raster in the same way to mimic that behavior
|
||||
pygeoprocessing.create_raster_from_vector_extents(
|
||||
base_proj_vector_path, initial_raster_path,
|
||||
(90, -90), wind_energy._TARGET_DATA_TYPE, wind_energy._TARGET_NODATA)
|
||||
# We have points at the four corners; update two of those values
|
||||
arr = pygeoprocessing.raster_to_numpy_array(initial_raster_path)
|
||||
arr[0, 0] = 1
|
||||
arr[-1, -1] = 1
|
||||
geometries = [
|
||||
Point(pos_x + 250, pos_y - 150), # out of extent
|
||||
Point(pos_x + 50, pos_y), # valid
|
||||
Point(pos_x + 150, pos_y), # invalid: nodata
|
||||
Point(pos_x + 250, pos_y), # out of extent
|
||||
Point(pos_x + 50, pos_y - 150), # invalid: nodata
|
||||
Point(pos_x + 150, pos_y - 150)] # valid
|
||||
|
||||
raster_info = pygeoprocessing.get_raster_info(initial_raster_path)
|
||||
origin = (raster_info['geotransform'][0], raster_info['geotransform'][3])
|
||||
base_vector_path = os.path.join(self.workspace_dir, 'base_vector.shp')
|
||||
pygeoprocessing.shapely_geometry_to_vector(
|
||||
geometries, base_vector_path, projection_wkt, 'ESRI Shapefile',
|
||||
fields=fields, attribute_list=attrs, ogr_geom_type=ogr.wkbPoint)
|
||||
|
||||
# Write modified raster values to use for pixel-picking
|
||||
final_raster_path = os.path.join(
|
||||
self.workspace_dir, 'modified_raster.tif')
|
||||
# Setup parameters for create raster
|
||||
matrix = numpy.array([
|
||||
[1, -1],
|
||||
[-1, 1]], dtype=numpy.int32)
|
||||
# Create raster to use for testing input
|
||||
raster_path = os.path.join(
|
||||
self.workspace_dir, 'raster.tif')
|
||||
pygeoprocessing.numpy_array_to_raster(
|
||||
arr, wind_energy._TARGET_NODATA, raster_info['pixel_size'], origin,
|
||||
raster_info['projection_wkt'], final_raster_path
|
||||
)
|
||||
matrix, -1, (100, -100), origin, projection_wkt,
|
||||
raster_path)
|
||||
|
||||
output_vector_path = os.path.join(
|
||||
self.workspace_dir, 'output_vector.shp')
|
||||
target_vector_path = os.path.join(self.workspace_dir, 'target_vector.shp')
|
||||
wind_energy._index_raster_values_to_point_vector(
|
||||
base_proj_vector_path, [(final_raster_path, 'TestVal')],
|
||||
output_vector_path, mask_keys=['TestVal'],
|
||||
mask_field="Masked"
|
||||
)
|
||||
base_vector_path, [(raster_path, 'TestVal')],
|
||||
target_vector_path, mask_keys=['TestVal'],
|
||||
mask_field="Masked")
|
||||
|
||||
# Confirm that masked-out points still exist in the base vector,
|
||||
# but have "Masked" values of 1
|
||||
base_vector = gdal.OpenEx(base_proj_vector_path, gdal.OF_VECTOR)
|
||||
# but have "Masked" values of 1. Out-of-extent points should
|
||||
# have values of None
|
||||
base_vector = gdal.OpenEx(base_vector_path, gdal.OF_VECTOR)
|
||||
base_layer = base_vector.GetLayer()
|
||||
self.assertEqual(base_layer.GetFeatureCount(), 4)
|
||||
id_masked = [(feat.GetField('pt_id'), feat.GetField('Masked'))
|
||||
self.assertEqual(base_layer.GetFeatureCount(), 6)
|
||||
id_masked = [(feat.GetField('id'), feat.GetField('Masked'))
|
||||
for feat in base_layer]
|
||||
self.assertEqual(id_masked, [(0, 0), (1, 1), (2, 1), (3, 0)])
|
||||
self.assertEqual(id_masked,
|
||||
[(0, None), (1, 0), (2, 1), (3, None), (4, 1), (5, 0)])
|
||||
base_layer = None
|
||||
base_vector = None
|
||||
|
||||
# Confirm that the target vector has only valid points
|
||||
target_vector = gdal.OpenEx(output_vector_path, gdal.OF_VECTOR)
|
||||
target_vector = gdal.OpenEx(target_vector_path, gdal.OF_VECTOR)
|
||||
target_layer = target_vector.GetLayer()
|
||||
self.assertEqual(target_layer.GetFeatureCount(), 2)
|
||||
feat_ids = [feat.GetField('pt_id') for feat in target_layer]
|
||||
self.assertEqual(feat_ids, [0, 3])
|
||||
feat_ids = [feat.GetField('id') for feat in target_layer]
|
||||
self.assertEqual(feat_ids, valid_ids)
|
||||
target_layer = None
|
||||
target_vector = None
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue