Carson Farmer
2015-10-22 19:13:54 UTC
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()
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()