Skip to content

Commit

Permalink
Merge branch 'main' of github.com:natcap/pygeoprocessing into feature…
Browse files Browse the repository at this point in the history
…/373-add-lrucache-to-convolve2d

Conflicts:
	HISTORY.rst
  • Loading branch information
phargogh committed Oct 11, 2024
2 parents 93993bb + ae34e3d commit 74a8086
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 15 deletions.
2 changes: 2 additions & 0 deletions HISTORY.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@ Unreleased Changes
------------------
* Removing the ``numpy<2`` constraint for requirements.txt that should have
been included in the 2.4.5 release. https://github.com/natcap/pygeoprocessing/issues/396
* Handling GDAL-based ``RuntimeError``s raised during ``pygeoprocessing.reproject_vector``.
https://github.com/natcap/pygeoprocessing/issues/409
* An LRU Cache has been added to ``convolve_2d`` which should reduce runtimes
on large rasters by approximately 10%. https://github.com/natcap/pygeoprocessing/issues/373

Expand Down
36 changes: 22 additions & 14 deletions src/pygeoprocessing/geoprocessing.py
Original file line number Diff line number Diff line change
Expand Up @@ -2243,21 +2243,29 @@ def reproject_vector(
100.0 * float(feature_index+1) / (layer.GetFeatureCount()),
os.path.basename(target_path))

geom = base_feature.GetGeometryRef()
if geom is None:
# we encountered this error occasionally when transforming clipped
# global polygons. Not clear what is happening but perhaps a
# feature was retained that otherwise wouldn't have been included
# in the clip
error_count += 1
continue
try:
geom = base_feature.GetGeometryRef()
if geom is None:
# we encountered this error occasionally when transforming
# clipped global polygons. Not clear what is happening but
# perhaps a feature was retained that otherwise wouldn't have
# been included in the clip
error_count += 1
continue

# Transform geometry into format desired for the new projection
error_code = geom.Transform(coord_trans)
if error_code != 0: # error
# this could be caused by an out of range transformation
# whatever the case, don't put the transformed poly into the
# output set
# Transform geometry into format desired for the new projection
error_code = geom.Transform(coord_trans)
if error_code != 0: # error
# this could be caused by an out of range transformation
# whatever the case, don't put the transformed poly into the
# output set
error_count += 1
continue
except RuntimeError as error:
# RuntimeError: GDAL's base error when geometries cannot be
# returned or transformed.
LOGGER.debug("Skipping feature %s due to %s",
feature_index, str(error))
error_count += 1
continue

Expand Down
44 changes: 43 additions & 1 deletion tests/test_geoprocessing.py
Original file line number Diff line number Diff line change
Expand Up @@ -498,7 +498,6 @@ def test_reproject_vector(self):
layer = None
vector = None


def test_reproject_vector_partial_fields(self):
"""PGP.geoprocessing: reproject vector with partial field copy."""
# Create polygon shapefile to reproject
Expand Down Expand Up @@ -798,6 +797,49 @@ def test_reproject_vector_utm_to_utm(self):
osr.SpatialReference(result_reference.ExportToWkt()).IsSame(
osr.SpatialReference(target_reference.ExportToWkt())))

def test_reproject_vector_handles_bad_data(self):
"""PGP.geoprocessing: reproject vector with bad data."""
vector_path = os.path.join(self.workspace_dir, 'bad_vector.shp')
driver = ogr.GetDriverByName('ESRI Shapefile')
vector = driver.CreateDataSource(vector_path)
srs = osr.SpatialReference()
srs.ImportFromEPSG(26710) # NAD27 / UTM zone 10N
layer = vector.CreateLayer(
'bad_layer', srs=srs, geom_type=ogr.wkbPoint)

# No/empty geometry
feature = ogr.Feature(layer.GetLayerDefn())
layer.CreateFeature(feature)

# Create a point at the centroid of the UTM zone 10N bounding box
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetGeometry(
ogr.CreateGeometryFromWkt('POINT (2074757.82 7209331.79)'))

layer = None
vector = None

# Our target UTM zone is 59S (in NZ), so the points should not be
# usable.
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(2134) # NZGD2000 / UTM zone 59S
target_srs_wkt = target_srs.ExportToWkt()

target_path = os.path.join(self.workspace_dir, 'target_vector.shp')
pygeoprocessing.reproject_vector(
base_vector_path=vector_path,
target_projection_wkt=target_srs_wkt,
target_path=target_path,
)

# verify that both fields were skipped.
try:
target_vector = gdal.OpenEx(target_path)
target_layer = target_vector.GetLayer()
self.assertEqual(target_layer.GetFeatureCount(), 0)
finally:
target_vector = target_layer = None

def test_calculate_disjoint_polygon_set(self):
"""PGP.geoprocessing: test calc_disjoint_poly no/intersection."""
gpkg_driver = ogr.GetDriverByName('GPKG')
Expand Down

0 comments on commit 74a8086

Please sign in to comment.