Skip to content

Commit

Permalink
modified transform geom to accept arrays of geoms
Browse files Browse the repository at this point in the history
  • Loading branch information
@snowman2
snowman2 committed Oct 30, 2019
1 parent d0e94b5 commit 1ae3dcf
Show file tree
Hide file tree
Showing 3 changed files with 131 additions and 64 deletions.
140 changes: 78 additions & 62 deletions fiona/_transform.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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(
<const OGRGeometry *>src_ogr_geom,
<OGRCoordinateTransformation *>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(
<const OGRGeometry *>src_ogr_geom,
<OGRCoordinateTransformation *>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)
Expand All @@ -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]
Expand All @@ -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]
Expand All @@ -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 = []
Expand All @@ -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
5 changes: 3 additions & 2 deletions fiona/transform.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`.
Expand All @@ -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.
Expand Down
50 changes: 50 additions & 0 deletions tests/test_transform.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 1ae3dcf

Please sign in to comment.