[Community] Tolerances

[Carson Farmer]

Hi all, I’m wondering if someone can provide some guidance on a ‘problem’
I’m working on. In the following example, I’m trying to ‘stitch’ together
geometries which cross the anti-meridian. Here I’m working with a (subset
of) the Russian geometries that comes from the Natural Earth 110m Admin
Shapefile (I’ve just included the WKT to make things more reproducible).
See this Gist <https://gist.github.com/carsonfarmer/e63d21dd3cf9105505db>
for the full example. The question is: “is it possible to use ‘fuzzy’
tolerances in any of the Shapely/GEOS geometry functions (or linemerge in
particular )?” or “is it possible to work around the fact that I can’t use
fuzzy tolerances?”. I don’t mind making copies of things or passing things
around to get this to work
 I’d just like to avoid having to decrease the
output precision of my merged geometries if possible

Alternatively, is there a better way to stitch geometries which cross the
anti-meridian using shapely? The TopoJSON commadline utility does it, but I
want to stay in Python land...

from shapely import wit
from shapely import ops
import matplotlib.pyplot as plt
%matplotlib inline

# These are provided in full in the Gist
a = wkt.loads("LINESTRING (-180 68.96363636363637

b = wkt.loads("LINESTRING (180 64.974329999999999


def plot_line(ob, ax, **kwargs):
x, y = ob.xy
ax.plot(x, y, linewidth=1, zorder=1, **kwargs)

def simplify(x, y, tol=8):
return round(x, tol), round(y, too)

def shift(x, y):
"""Sometimes x is a tuple, sometimes a float?!"""
try:
x = [i + 360 if i < 0 else i for i in x]
except:
x = x + 360 if x < 0 else x
return x, y

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4))

# This obviously isn't going to work

l = ops.linemerge([a, b])
for i in l:
plot_line(i, ax1)

# This should work, but the points aren't *exactly* coincident
l = ops.linemerge([ops.transform(shift, a), b])
for i in l:
plot_line(i, ax2)

# This *does* work, but the process is ‘destructive'
# (i.e., we lose precision)
# You can check this by looking at l.coords
l = ops.linemerge([ops.transform(simplify, ops.transform(shift, i)) for i in
[a, b]])
plot_line(l, ax3)

plt.show()

[Sean Gillies]

Hi Carson,

I've scoped out the work of adding fixed precision geometry models to
Shapely but haven't begun. It's a pretty big lift and I've got too
much other stuff on my plate right now. For now, Shapely has only a
floating point precision model. This causes many problems with Natural
Earth: that dataset is full of features that are invalid under the
floating point model (due to microscopic self-intersections &c) but
that would be perfectly valid with a fixed precision model.

For merging anti-meridian crossers: once you've shifted geometry parts
from one side by 360 degrees, you should be able to use
shapely.ops.unary_union.

Hi all, I’m wondering if someone can provide some guidance on a ‘problem’
I’m working on. In the following example, I’m trying to ‘stitch’ together
geometries which cross the anti-meridian. Here I’m working with a (subset
of) the Russian geometries that comes from the Natural Earth 110m Admin
Shapefile (I’ve just included the WKT to make things more reproducible). See
this Gist for the full example. The question is: “is it possible to use
‘fuzzy’ tolerances in any of the Shapely/GEOS geometry functions (or
linemerge in particular )?” or “is it possible to work around the fact that
I can’t use fuzzy tolerances?”. I don’t mind making copies of things or
passing things around to get this to work… I’d just like to avoid having to
decrease the output precision of my merged geometries if possible…
Alternatively, is there a better way to stitch geometries which cross the
anti-meridian using shapely? The TopoJSON commadline utility does it, but I
want to stay in Python land...
from shapely import wit
from shapely import ops
import matplotlib.pyplot as plt
%matplotlib inline
# These are provided in full in the Gist
a = wkt.loads("LINESTRING (-180 68.96363636363637…
b = wkt.loads("LINESTRING (180 64.974329999999999…
x, y = ob.xy
ax.plot(x, y, linewidth=1, zorder=1, **kwargs)
return round(x, tol), round(y, too)
"""Sometimes x is a tuple, sometimes a float?!"""
x = [i + 360 if i < 0 else i for i in x]
x = x + 360 if x < 0 else x
return x, y
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4))
# This obviously isn't going to work…
l = ops.linemerge([a, b])
plot_line(i, ax1)
# This should work, but the points aren't exactly coincident
l = ops.linemerge([ops.transform(shift, a), b])
plot_line(i, ax2)
# This does work, but the process is ‘destructive'
# (i.e., we lose precision)
# You can check this by looking at l.coords
l = ops.linemerge([ops.transform(simplify, ops.transform(shift, i)) for i in
[a, b]])
plot_line(l, ax3)
plt.show()
_______________________________________________
Community mailing list
http://lists.gispython.org/mailman/listinfo/community

[Carson Farmer]

Hi Sean, thanks for the response. This is pretty much what I figured :p
I'll just stick with the 'fixed precision like hack' and use that as a
workaround. It's not the end of the world (get it?)... I guess I could
store a mapping of the floating point to the rounded values (not ideal) if
I really needed to get them back...

Cheers,
Carson

Post by Sean Gillies
Hi Carson,
I've scoped out the work of adding fixed precision geometry models to
Shapely but haven't begun. It's a pretty big lift and I've got too
much other stuff on my plate right now. For now, Shapely has only a
floating point precision model. This causes many problems with Natural
Earth: that dataset is full of features that are invalid under the
floating point model (due to microscopic self-intersections &c) but
that would be perfectly valid with a fixed precision model.
For merging anti-meridian crossers: once you've shifted geometry parts
from one side by 360 degrees, you should be able to use
shapely.ops.unary_union.

Post by Carson Farmer
Hi all, I’m wondering if someone can provide some guidance on a ‘problem’
I’m working on. In the following example, I’m trying to ‘stitch’ together
geometries which cross the anti-meridian. Here I’m working with a (subset
of) the Russian geometries that comes from the Natural Earth 110m Admin
Shapefile (I’ve just included the WKT to make things more reproducible).

See

Post by Carson Farmer
this Gist for the full example. The question is: “is it possible to use
‘fuzzy’ tolerances in any of the Shapely/GEOS geometry functions (or
linemerge in particular )?” or “is it possible to work around the fact

that

Post by Carson Farmer
I can’t use fuzzy tolerances?”. I don’t mind making copies of things or
passing things around to get this to work
 I’d just like to avoid having

to

Post by Carson Farmer
decrease the output precision of my merged geometries if possible

Alternatively, is there a better way to stitch geometries which cross the
anti-meridian using shapely? The TopoJSON commadline utility does it,

but I

Post by Carson Farmer
want to stay in Python land...
from shapely import wit
from shapely import ops
import matplotlib.pyplot as plt
%matplotlib inline
# These are provided in full in the Gist
a = wkt.loads("LINESTRING (-180 68.96363636363637

b = wkt.loads("LINESTRING (180 64.974329999999999

x, y = ob.xy
ax.plot(x, y, linewidth=1, zorder=1, **kwargs)
return round(x, tol), round(y, too)
"""Sometimes x is a tuple, sometimes a float?!"""
x = [i + 360 if i < 0 else i for i in x]
x = x + 360 if x < 0 else x
return x, y
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4))
# This obviously isn't going to work

l = ops.linemerge([a, b])
plot_line(i, ax1)
# This should work, but the points aren't exactly coincident
l = ops.linemerge([ops.transform(shift, a), b])
plot_line(i, ax2)
# This does work, but the process is ‘destructive'
# (i.e., we lose precision)
# You can check this by looking at l.coords
l = ops.linemerge([ops.transform(simplify, ops.transform(shift, i)) for

i in

Post by Carson Farmer
[a, b]])
plot_line(l, ax3)
plt.show()
_______________________________________________
Community mailing list
http://lists.gispython.org/mailman/listinfo/community

--
Sean Gillies
_______________________________________________
Community mailing list
http://lists.gispython.org/mailman/listinfo/community