Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

modified transform geom to accept arrays of geoms #811

Merged
merged 2 commits into from
Nov 5, 2019
Merged

modified transform geom to accept arrays of geoms #811

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a11d4a9
modified transform geom to accept arrays of geoms
snowman2 Oct 12, 2019
3da7b85
Merge branch 'maint-1.9' into transform_array
sgillies Nov 5, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
140 changes: 78 additions & 62 deletions fiona/_transform.pyx
View file
Open in desktop
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):
snowman2 marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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 (or a list of GeoJSON-like geometries
if an iterable was given as input) with transformed coordinates. 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
77 changes: 45 additions & 32 deletions tests/test_transform.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,49 +9,62 @@
from .conftest import requires_gdal_lt_3


@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": [
TEST_GEOMS = [
{"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],
]
],
},
{
"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],
]
]
],
},
],
)
]
],
},
]


@pytest.mark.parametrize("geom", TEST_GEOMS)
def test_transform_geom_with_z(geom):
"""Transforming a geom with Z succeeds"""
transform.transform_geom("epsg:4326", "epsg:3857", geom, precision=3)


@pytest.mark.parametrize("geom", TEST_GEOMS)
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


@requires_gdal_lt_3
def test_transform_geom_null_dest():
failed_geom = {
Expand Down