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()
_______________________________________________
Community mailing list
Community@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community
