From 1ae3dcf12c3efd09bd9d4e67c76f3d0e2951dc3b Mon Sep 17 00:00:00 2001 From: snowman2 Date: Fri, 11 Oct 2019 19:50:00 -0500 Subject: [PATCH] modified transform geom to accept arrays of geoms --- fiona/_transform.pyx | 140 ++++++++++++++++++++++------------------ fiona/transform.py | 5 +- tests/test_transform.py | 50 ++++++++++++++ 3 files changed, 131 insertions(+), 64 deletions(-) diff --git a/fiona/_transform.pyx b/fiona/_transform.pyx index 61ee42450..0ba48c3c1 100644 --- a/fiona/_transform.pyx +++ b/fiona/_transform.pyx @@ -11,7 +11,7 @@ from fiona cimport _cpl, _crs, _csl, _geometry from fiona._crs cimport OGRSpatialReferenceH from fiona._shim cimport osr_set_traditional_axis_mapping_strategy -from fiona.compat import UserDict +from fiona.compat import UserDict, DICT_TYPES cdef extern from "ogr_geometry.h" nogil: @@ -115,63 +115,38 @@ def _transform(src_crs, dst_crs, xs, ys): return res_xs, res_ys -def _transform_geom( - src_crs, dst_crs, geom, antimeridian_cutting, antimeridian_offset, - precision): - """Return a transformed geometry.""" - cdef char *proj_c = NULL - cdef char *key_c = NULL - cdef char *val_c = NULL - cdef char **options = NULL - cdef OGRSpatialReferenceH src = NULL - cdef OGRSpatialReferenceH dst = NULL - cdef void *transform = NULL - cdef OGRGeometryFactory *factory = NULL +cdef object _transform_single_geom( + object single_geom, + OGRGeometryFactory *factory, + void *transform, + char **options, + object precision +): cdef void *src_ogr_geom = NULL cdef void *dst_ogr_geom = NULL cdef int i - if src_crs and dst_crs: - src = _crs_from_crs(src_crs) - dst = _crs_from_crs(dst_crs) - transform = _crs.OCTNewCoordinateTransformation(src, dst) - # Transform options. - options = _csl.CSLSetNameValue( - options, "DATELINEOFFSET", - str(antimeridian_offset).encode('utf-8')) - if antimeridian_cutting: - options = _csl.CSLSetNameValue(options, "WRAPDATELINE", "YES") - - factory = new OGRGeometryFactory() - src_ogr_geom = _geometry.OGRGeomBuilder().build(geom) - dst_ogr_geom = factory.transformWithOptions( - src_ogr_geom, - transform, - options) - if dst_ogr_geom == NULL: - g = None - warnings.warn( - "Full reprojection failed, but partial is possible. To enable partial " - "reprojection wrap the transform_geom call like so:\n" - "with fiona.Env(OGR_ENABLE_PARTIAL_REPROJECTION=True):\n" - " transform_geom(...)" - ) - else: - g = _geometry.GeomBuilder().build(dst_ogr_geom) - _geometry.OGR_G_DestroyGeometry(dst_ogr_geom) - if src_ogr_geom != NULL: - _geometry.OGR_G_DestroyGeometry(src_ogr_geom) - _crs.OCTDestroyCoordinateTransformation(transform) - if options != NULL: - _csl.CSLDestroy(options) - _crs.OSRRelease(src) - _crs.OSRRelease(dst) - + src_ogr_geom = _geometry.OGRGeomBuilder().build(single_geom) + dst_ogr_geom = factory.transformWithOptions( + src_ogr_geom, + transform, + options) + if dst_ogr_geom == NULL: + out_geom = None + warnings.warn( + "Full reprojection failed, but partial is possible. To enable partial " + "reprojection wrap the transform_geom call like so:\n" + "with fiona.Env(OGR_ENABLE_PARTIAL_REPROJECTION=True):\n" + " transform_geom(...)" + ) else: - g = geom - - if precision >= 0 and g is not None: - if g['type'] == 'Point': - coords = list(g['coordinates']) + out_geom = _geometry.GeomBuilder().build(dst_ogr_geom) + _geometry.OGR_G_DestroyGeometry(dst_ogr_geom) + if src_ogr_geom != NULL: + _geometry.OGR_G_DestroyGeometry(src_ogr_geom) + + if out_geom is not None and precision >= 0: + if out_geom['type'] == 'Point': + coords = list(out_geom['coordinates']) x, y = coords[:2] x = round(x, precision) y = round(y, precision) @@ -180,8 +155,8 @@ def _transform_geom( z = coords[2] new_coords.append(round(z, precision)) - elif g['type'] in ['LineString', 'MultiPoint']: - coords = list(zip(*g['coordinates'])) + elif out_geom['type'] in ['LineString', 'MultiPoint']: + coords = list(zip(*out_geom['coordinates'])) xp, yp = coords[:2] xp = [round(v, precision) for v in xp] yp = [round(v, precision) for v in yp] @@ -192,9 +167,9 @@ def _transform_geom( else: new_coords = list(zip(xp, yp)) - elif g['type'] in ['Polygon', 'MultiLineString']: + elif out_geom['type'] in ['Polygon', 'MultiLineString']: new_coords = [] - for piece in g['coordinates']: + for piece in out_geom['coordinates']: coords = list(zip(*piece)) xp, yp = coords[:2] xp = [round(v, precision) for v in xp] @@ -206,8 +181,8 @@ def _transform_geom( else: new_coords.append(list(zip(xp, yp))) - elif g['type'] == 'MultiPolygon': - parts = g['coordinates'] + elif out_geom['type'] == 'MultiPolygon': + parts = out_geom['coordinates'] new_coords = [] for part in parts: inner_coords = [] @@ -224,6 +199,47 @@ def _transform_geom( inner_coords.append(list(zip(xp, yp))) new_coords.append(inner_coords) - g['coordinates'] = new_coords + out_geom['coordinates'] = new_coords + return out_geom + + +def _transform_geom( + src_crs, dst_crs, geom, antimeridian_cutting, antimeridian_offset, + precision): + """Return a transformed geometry.""" + cdef char *proj_c = NULL + cdef char *key_c = NULL + cdef char *val_c = NULL + cdef char **options = NULL + cdef OGRSpatialReferenceH src = NULL + cdef OGRSpatialReferenceH dst = NULL + cdef void *transform = NULL + cdef OGRGeometryFactory *factory = NULL + if not all([src_crs, dst_crs]): + raise RuntimeError("Must provide a source and destination CRS.") + src = _crs_from_crs(src_crs) + dst = _crs_from_crs(dst_crs) + transform = _crs.OCTNewCoordinateTransformation(src, dst) + # Transform options. + options = _csl.CSLSetNameValue( + options, "DATELINEOFFSET", + str(antimeridian_offset).encode('utf-8')) + if antimeridian_cutting: + options = _csl.CSLSetNameValue(options, "WRAPDATELINE", "YES") + + factory = new OGRGeometryFactory() - return g + if isinstance(geom, DICT_TYPES): + out_geom = _transform_single_geom(geom, factory, transform, options, precision) + else: + out_geom = [ + _transform_single_geom(single_geom, factory, transform, options, precision) + for single_geom in geom + ] + + _crs.OCTDestroyCoordinateTransformation(transform) + if options != NULL: + _csl.CSLDestroy(options) + _crs.OSRRelease(src) + _crs.OSRRelease(dst) + return out_geom \ No newline at end of file diff --git a/fiona/transform.py b/fiona/transform.py index b380b6aff..552021ef8 100644 --- a/fiona/transform.py +++ b/fiona/transform.py @@ -58,7 +58,7 @@ def transform_geom( on the "destination" or "to" side of the transformation. geom: obj A GeoJSON-like geometry object with 'type' and 'coordinates' - members. + members or an iterable of GeoJSON-like geometry objects. antimeridian_cutting: bool, optional ``True`` to cut output geometries in two at the antimeridian, the default is ``False`. @@ -72,7 +72,8 @@ def transform_geom( Returns ------- obj - A new GeoJSON-like geometry with transformed coordinates. Note + A new GeoJSON-like geometry with transformed coordinates or a list + of GeoJSON-like geometries if a list was given as input. Note that if the output is at the antimeridian, it may be cut and of a different geometry ``type`` than the input, e.g., a polygon input may result in multi-polygon output. diff --git a/tests/test_transform.py b/tests/test_transform.py index 89addc99a..0ac194ffd 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -74,3 +74,53 @@ def test_transform_geom_null_dest(): antimeridian_cutting=True, precision=2, ) is None + + +@pytest.mark.parametrize( + "geom", + [ + {"type": "Point", "coordinates": [0.0, 0.0, 1000.0]}, + { + "type": "LineString", + "coordinates": [[0.0, 0.0, 1000.0], [0.1, 0.1, -1000.0]], + }, + { + "type": "MultiPoint", + "coordinates": [[0.0, 0.0, 1000.0], [0.1, 0.1, -1000.0]], + }, + { + "type": "Polygon", + "coordinates": [ + [ + [0.0, 0.0, 1000.0], + [0.1, 0.1, -1000.0], + [0.1, -0.1, math.pi], + [0.0, 0.0, 1000.0], + ] + ], + }, + { + "type": "MultiPolygon", + "coordinates": [ + [ + [ + [0.0, 0.0, 1000.0], + [0.1, 0.1, -1000.0], + [0.1, -0.1, math.pi], + [0.0, 0.0, 1000.0], + ] + ] + ], + }, + ], +) +def test_transform_geom_array_z(geom): + """Transforming a geom array with Z succeeds""" + g2 = transform.transform_geom( + "epsg:4326", + "epsg:3857", + [geom for _ in range(5)], + precision=3, + ) + assert isinstance(g2, list) + assert len(g2) == 5